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ABSTRACT 


This work describes the effects of erosion on the heat transfer characteristics on thrust vector 
control vanes exposed to aluminized propellant exhaust flows. This was accomplished using an 
inverse heat transfer parameter identification of quarter scale models. The model is based on a four 
node lumped parameter system with two heat energy inputs. The erosion is modeled as decreasing 
the geometric dimensions linearly as a function of time and the percentage of aluminum in the 
propellant. Excellent agreement was found between experimental and model temperature profiles. 


The heat transfer coefficients of the vanes were found to decrease with increasing erosion rates. 
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I. INTRODUCTION 


This thesis is a continuation of work done by the Naval 
Air Warfare Center Weapons Division (NAWCWPNS) and thesis work 
at the Naval Postgraduate School (NPS) to provide a better 
understanding of the heat transfer characteristics of jet 
vanes used for thrust vector control (TVC) of vertical launch 
missiles. This is accomplished using an inverse heat transfer 
parameter identification of quarter scale replicas which can 
be used to find full scale results. 

Thrust vector control is a process by which jet vanes are. 
inserted into the exhaust plume of a missile to control the 
flight path prior to the missile obtaining the required 
velocity for the external control surfaces to take effect. 


{[Ref. 1] A schematic for the TVC system is shown in Figure 1. 
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Figure 1 Thrust Vector Control System Schematic 


Due to the harsh thermal environment that the vanes are 
exposed to, a better understanding of the heat transfer 
processes which take place will help in the improved design of 
jet vanes. This will lead to longer operation and the ability 
to use propellants that burn hotter and use a higher 
percentage of aluminum for greater momentum flux and better 
performance. [Ref. 2:p. 1] 

There are five basic steps in determining the heat 
transfer characteristics of the vane: 

A Develop a mathematical model of the heat transfer 
processes which take place in the vane. It is expressed in 
terms of a number of physical constants, some of which are 
known, some of which are to be determined. [Ref. 3:p. 1] 

2. Gather experimental data in the form of temperature- 
time data at selected locations on the vane. 

3. Compare the predicted and experimental temperature-time 
data. 

4. Use the differences between the simulated and actual 
temperatures to drive a systematic adjustment of unknown model 
parameters in an optimization routine. The process is 
repeated until the experimental and theoretical data 
differences are minimized in a least-squares sense. 

{Ref. 3:p. 2] 

5. Calculate the heat transfer parameters of the system 

using the physical parameters of the model which give the best 


estimate of the actual behavior. 


Previous work has concentrated on using parametric system 
identification to validate the use of full and quarter scale 
models to predict the heat transfer characteristics for full 
scale vanes in a non-erosive environment. The research in 
this report extends the quarter scale model to an erosive 


environment. 


It. THEORY 


A. BACKGOUND 
1. Physical Description 

The main pieces of equipment used in the experimental 
tests are the rocket motor and the jet vanes. The rocket 
motor is set up to provide a constant thrust-time profile. 
The propellants used in the motor are aluminized hydroxyl- 
terminated polybutadiene (HTPB) with either 0%, 9%, or 18% Al 
by weight. The jet vanes are made from pressed and sintered 
tungstan powder that is infiltrated by 10% copper by weight.. 
There are four vanes for each motor. The experimental setup 


is shown in Figure 2. [Ref. 2:p. 1,2] 
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Figure 2 Vane and Motor Assembly for Experimental Tests 


The experimental tests are conducted as either full or 
quarter scale. The quarter scale tests have several 
advantages. Most important is the cost savings over a full 
scale test. The reduced size of the motor, vane and test 
equipment account for much of the savings. [Ref. 4:p. 15,16] 
The biggest disadvantage of the quarter scale vane comes in 
the placement of the thermocouples. Whereas in the full scale 
vane the thermocouples can be placed inside the vane, for the 
quarter scale vane the thermocouples must be placed on the 


vane shaft. The thermocouple placement is contrasted in 


Figure 3. 


Figure 3 Thermocouple Placement for Full and Quarter Scale 


2. Basic Modeling 
“In order to predict the thermal response of the jet 


vane, a simple model had to be developed. The model had to 


consider the physical characteristics of the vane and the heat 
transfer processes that were taking place. 

The physical quantities can be broken into two 
categories: material and geometric. The material properties 
considered were the vane density, thermal conductivity, and 
specific heat. The geometric properties considered were the 
conductive lengths, cross sectional and surface areas and 
volume. 

The heat transfer processes considered were 
convection at the surface of the vane and conduction of heat 


through the vane. 


Figure 4 Thermal Energy Node Model 


“Heat transfer in the vane is modeled by applying the 
law of conservation of energy. Energy balance equations can 


be derived using a model consisting of thermal resistances and 


capacitances driven by the temperature difference between the 


The energy balance for Figure 4 is, 


nodes. 
(1) 


T,= - 
* GR; CR; OR, CR, 


The convective resistance is found by, 


where Tj>T,>Ty- 


1 
(2) 


nA, 


R= 


where h is the convective heat transfer coefficient and A, is 


the surface area. The conductive resistance is found by, 
£L 
(3) 


7 
S 


the conductive length, k is the thermal 


where L is 
The thermal 


conductivity and A, is the cross sectional area. 


capacitance is given by, 
(4) 


C=pVC, 


where p is the material density, V is the volume and C, is the 


material specific heat. 


Lumped Parameter 
The nodes of the basic model lend itself to dividing 
For the full 


3. 


the vane into different sections, or lumps. 


scale model, the vane was geometrically divided into three 
separate sections: the tip, fin and shaft. A node is located 
at the center of each section. The sections are defined as 
shown in Figure 5. For the quarter scale model, a fourth node 
was added at the mount to account for the different 
thermocouple placement. 
4. PSI Process 

A simple model was needed that could easily be changed 
for different materials, geometries and exhaust conditions. 
This lead to parametric system identification (PSI). PSI is 
a computer based procedure where the parameters of a model are 
changed until a best fit approximation in a least squares 
sense to experimental data is obtained. [Ref. 5:p. 6] 

Parameter identification has several advantages over 
the other modeling choice, computational fluid dynamics (CFD). 
Creating a mathmatical model of the vane using CFD is almost 
impossible due to the complexity of the exhaust flow. The jet 
vane must operate in a high temperature, three-dimensional, 
turbulent, compressible supersonic flow. [Ref. 5:p. 3,4] PSI 
ignores these complexities and focuses on the end result. 
This makes PSI not only simpler, but the information that 
comes out of the PSI model can easily be used in improving the 


design. PSI also handles nonlinear conditions such as 


ablation. [Ref. 6: p.2] 


Figure 5 Jet Vane as a Lumped Parameter Model 
(Actual Design Indicated by Dotted Lines) 


The simple three node model shown in Figure 6 can 
serve as baseline for other models. The model is driven by 
two heat sources, represented by temperatures T,, and T,,, which 
are the stagnation and free stream recovery temperatures, 
respectively. The temperature T,, heats the vane through 


convective heat transfer at node one through the thermal 


Figure 6 Three Node Jet Vane Model 


resistance R,, and stores the energy as a thermal capacitance 
an, iC). The same process occurs at node two with recovery 
temperature T,,, thermal resistance R,, and thermal capacitance 
C,. Node three stores energy in thermal capacitance C, and is 
connected to ground through thermal resistance R,. All nodes 
are coupled by conductive resistances. [Ref. 6:p. 4] Applying 
the law-of conservation of energy to the system leads to the 


following equations: 


10 


fs Tt 7; oe T; 
. CR, CyRo3 CRyg 
letting, 
1 1 
Q,5= a,,= 
ne CLRy2 es CoRy2 
= 1 a.,= 1 
2° CaRos = C3R23 
1 
a b,,= 
7 C3R3¢ o Cy Rey 
1 
b= 
C2 Rpz 


Combining coefficients at the same temperatures gives, 


Ay, = 4,2 +B, 
227421 +42,+b,, 


433 74s2 Tyg 


Rewriting the equations, 


Ty =" 4,17, +4,2T, +B, ,T py 
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(5) 


(6) 


(7) 


(8) 


(9) 


(10) 


(11) 


(12) 


(13) 


(14) 


(15) 


T= G17, - 4327, +247; + bz 2T pz (16) 


T,=a,2T,-4,,T, (17) 


Rewriting into state-space form, T = AT + Bu, or 


P 7, 
kg Tay, Gar 0 7 Bb. 9 0 Tat 
T, = A, "422 423 T, + 0 bo 0 Tre (18) 
: -a 0 
T; Q 43a = T; r a 0 


The energy balance equations are a set of linear, ordinary 
differential equations which can be readily solved on a 
computer. This was done in a Fortran program using an IMSL 
subroutine called DIVPRK. DIVPRK solves a double precision 
initial value problem for ordinary differential equations 
using fifth-order and sixth-order Runge-Kutta-Verner methods. 
DIVPRK requires a user supplied subroutine called FCN which 
defines the set of equations to be solved. 

The main program containing DIVPRK and FCN is called 
SIM.FOR, and simulates the temperatures of the three node 
model. The model is driven by an input vector u which is the 
product of the recovery temperatures T,, and T, and a step 
function simulating the thrust. Physical and geometric data 
was used to calculate the internal thermal conductive 
resistances and capacitances which lead to coefficients in the 


A matrix. Since the inputs at nodes one and two from 


12 


convection and node three from ground are unknown, values for 
these resulting coefficients must be guessed. The output of 
the program is temperature-time data which is written to a 
data file called TEMP.MAT. This data can then be read into 
MATLAB and plotted. The purpose is to try to match calculated 
temperatures with known experimental temperature data at node 
two and validate the numerical approach. The results are 
shown in Figure 7. Although the node two temperatures are 
close, they are not identical. By extending the program to 
include an optimizer that could adjust the unknown A and B 
coefficients, a closer approximation could be found. 

This was done in a Fortran program called NODE3.FOR. 
It is in this parameter identification, or PID, program that 
the differential equations are set up and solved. First, 
physical and geometric data is read in froma data file called 
INPUT.DAT. This information is used to calculate the internal 
thermal conductive resistances and capactitances which lead to 
coefficients in the A matrix. Since the inputs at nodes one 
and two from convection and node three from ground are 
unknown, the resulting unknown coefficients from the A and B 
matrices are sent to the optimizer as variables to be found. 
The optimizer used is an IMSL routine called DBCLSF which uses 
a modified Levenberg-Marquardt method and an active set 
strategy to minimize an error in a least-squares sense subject 
to simple constraints placed on the variables by the user. 


DBCLSF calls a user written subroutine called TEMP that 


13 


THREE NODE MODEL _ 


Pod 
Me 
n 
2 
a 
. 

] 
a 
AJ 

~~ 
i] 
eZ 
2 
& 
t 
= 
sy 
Oy 
= 
eI 
be 


2.5 


TIME (sec) 


Figure 7 Simulated and Experimental Node Temperatures 


calculates the temperature-time history using the current 
parameters supplied by DBCLSF called from the PID program. It 
does this through DIVPRK and FCN. Once the temperature-time 
history -is calculated, an error function is returned to DBCLSF 
based on the differences between predicted and experimental 


temperature-time histories. The optimizer then adjusts the 
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unknown parameters and the process repeats until certain 


convergence criteria is met. 


B. PREVIOUS MODELS 

Work on the jet vane thermal model began at the Naval 
Postgraduate School (NPS) by Nunn and Kelleher [Ref. 7] in 
1986. Further development of the model was continued by Nunn 
(Ref. 5] and Hatzenbuehler. [Ref. 8] Hatzenbuehler was able 
to create a four node quarter scale model using PSI procedures 
and a computer software package called Matrix X. Reno [Ref. 
1] followed Hatzenbuehler and refined the four node model and 
attempted to compare the quarter scale results to full scale 
vanes, but was unsuccessful. More recent work has been done 
by Parker [Ref. 4]. He obtained good results using a full 
scale model of the jet vane. He also looked more closely at 
the scaling of the models and the applicability of quarter 
scale results to full scale vanes. He also found that 
existing quarter scale models did not provide an accurate 
picture of the heat transfer processes in the full scale 


vanes. 


C. THREE NODE FULL SCALE MODEL 
Parker’s five node full scale model was reduced to a three 
node full scale model to investigate whether the three fin 


nodes could be reduced to one node and obtain the same 
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results. Parker’s five node full scale model is shown in 


Figure 8. 


it | 
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Figure 8 Parker Five Node Full Scale Model 


The three node model was driven using the geometric data 
given in Table 1 and the following material data: p = 18310 
kg/m, k = 173 W/mK, and C, = 146 J/kgKk. The recovery 
temperatures used to drive the system were T,, = 2670 K and Try, 
= 2570 K. [Ref. 6:p. 7,8] These temperatures were contained 
in the input vector u, whose values were the product of the 
recovery temperature and a step function simulating the thrust 


function. 
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Table I GEOMETRIC DATA FOR FULL SCALE VANES 


The program found the values for b, = 1.0029 and b, = 
0.0809. This corresponds to the convection heat transfer 
coefficients of 16,025 W/mK and 1003 W/mK at the tip and fin 
respectively. The ground resistance was found to be 0.0001. 
These values were found to be reasonably close to those from 
Parker’s five node model. He found b,, = 1.3787, b, = 0.0862,. 
and the ground resistance to be 0.0001, while the convection 
heat transfer coefficients were 22027.5 W/mK and 1057 W/mK at 
the tip and fin respectively. [Ref. 4:p. 63] The temperature- 


time histories for both models are shown in Figure 9. 
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SIMULATION FOR THREE NODE FULL SCALE MODEL 
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Figure 9 Temperature-time Histories for Three and Five Node 
Full Scale Models 
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IIt. ABLATION EFFECTS 


A. ABLATION MODELING 

There was erosion in the quarter scale vanes exposed to 
aluminized propellant exhaust flows. For the 0% aluminized 
case, only 1% of the vane’s mass was lost. But for the 9% and 
18% aluminized cases, the loss became much more substantial. 
For the 9% case, 8% of the vane’s mass was lost. For the 18% 
aluminized case, 50% of the vane’s mass was lost. Vane mass 
loss was found to be nonlinear with the percentage of Al in 
the propellant. The relationship using an exponential 


function by an empirical fit was found to be 


mass loss=1.042e (9-2173) (tal) (19) 


Vane erosion profiles for the three cases are shown in Figure 
10. (Ref. 2:p. 6,7] At least part of this erosion was likely 
caused by ablation. Ablation is due to the melting of the 
surface of the vane [Ref. 9:p. 122]. 

A short FORTRAN program, COEF.FOR, was written to see how 
the known A matrix coefficients were affected by the mass 
loss. The geometric dimensions of length, area, and volume 
were modeled as decreasing linearly as a function of time and 
percent mass loss. The results for the 9% Al and 18% Al cases 


are shown in Figure 11. 
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Figure 10 Vane Erosion Profiles 


Several trends in Figure 11 are worth noting. The 
dominent coefficient in both cases is a,. This is expected 


since 


i: 
12> (20) 
2 GR 


and C, is small due to the small volume at the tip of the 
vane. Also note that the coefficients are nonlinear over time 


and the nonlinearity increases with increased mass loss. 


B. FOUR NODE QUARTER SCALE MODEL 
The erosion present in the quarter scale vanes when the 


aluminized propellant was used needed to be investigated. A 


four node quarter scale model had already been derived by 
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9% al CASE 
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Figure 11 Effect of Erosion on A Matrix Coefficients in the 
9% Al and 18% Al Cases 
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Reno. [Ref. 1] Application of the law of conservation of 


energy led to the following equations: 


T,=-a,,T,+4,,T,+b,,Tpy (21) 
T,=€2,T, ~ p27, *Az3T3+b22T a2 (22) 
T= 855 Tp ~243T, +4547, (23) 

Ty =4q3T,-2y4Ty (24) 


These equations needed to be modified though, since they 
did not include the effects of erosion. Erosion of the vane 
caused the geometric dimensions of the vane to change, while 
the material properties of density, thermal conductivity and 
specific heat remained constant. The program COEF.FOR modeled 
the changing geometric dimensions with time. All that was 
needed was to attach COEF.FOR to the main PID program as a 
subprogram. 

The other aspect of interest in the cases with aluminized 
propellant was whether the convective heat transfer 
coefficients were time varient. Once the values of b, and b, 
are found in the PID program, the program COEF.FOR can be 
modified so that the heat transfer coefficients can be 


calculated at every time step since 


- 1 ee 
Be RrAse ns Rp2Ase oe 
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1 1 
Ro;= R..= 
si bycy si B.C, ae 


and cC,, C,, A,, and A, are all time dependant. 


Cc. CONVERGING QUARTER SCALE MODEL WITH ABLATION 
1. Case 1: 0% Al in Propellant 

For case 1, data was taken for three seconds before 
thrust began to tailoff. This allowed for 61 temperature-time 
data points to be taken, or 20 per second. The data points on 
the vane corresponded to nodes three and four of the model. 
This data was read into the PID program NODE40.FOR along with 
the geometric data and the recovery temperatures. In the 
subroutine FCN, a delay of 0.3 seconds was used to account for 
the time before the thrust reached its steady state value. 
The results obtained were excellent; the square root of the 
sum of the squares of the difference between experimental and 
model temperatures at nodes three and four was only 1.19 
degrees Kelvin. A plot of the experimental and model 
temperatures is shown in Figure 12. 

The values obtained for the unknown variables were 
Aay=0.5376, ag=0.1528, ag=-0.1651, by=7.6511, and by,=0.0722. 
These variables led to resistance values of R,,=1.2221, 
Rpy=6.4795, and Rg=-1.8206. The negative value obtained for Ry 
indicates heating of the vane from ground. The convection 


heat transfer coefficients were calculated to be 30405.43 
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Figure 12 Case 1: Experimental and Model Temperatures Vs. 
Time 


W/mK and 222.43 W/mK at the tip and fin respectively. 
2. Case 2: 9% Al in Propellant 

The same procedure was done for case 2. Temperature- 
time data was only taken for two seconds before thrust 
tailoff. A delay of 0.7 seconds was used to account for the 
time before the thrust reached its steady state value. Again, 
the results were excellent: the sum of the squares difference 
was only 0.73 degrees Kelvin. A plot of the experimental and 
model temperatures is shown in Figure 13. 

‘The values obtained for the unknown variables were 


A@y=-0.2000, Aag=3.5088, A, =100.00, b,=3.4427, and b,=0.0837. 
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Figure 13 Case 2: Experimental and Model Temperatures Vs. 
Time 


The variables lead to resistance values of R,,=2.9135, 
Rp=5.9930, and R,=-0.1989. Again, the negative resistance of 
Rw indicates heating of the vane from ground. The convection 
heat transfer coefficients were calculated to be 13354.40 and 
251.80 at the tip and fin respectively. 

The values for b, and b, found from NODE49.FOR were 
added to the geometric and material data in COEF.FOR in order 
to calculate the convective heat transfer coefficients at 
every time step. The heat transfer coefficient at the tip 
decreased from an initial value of 13742.78 to the final value 
of 13354.40. The heat transfer coefficient for the fin 


decreased from an initial value of 259.17 to the final value 
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of 251.80. In both cases, there was only a three percent 
decrease. The coefficients are plotted versus time in ligure 


14. 


9% Al CASE 


Figure 14 Convective Heat Transfer Coefficients Plotted Vs. 
Time For Case 2. 


3. Case 3: 18% Al in Propellant 
The same procedure was done for case 3. Temperature- 
time data was only taken for 1.6 seconds before the severity 
of the erosion caused direct plume impingment to the vane 
shaft. (Ref. 2,p.9] A delay of 0.1 seconds was used to 
accent. fox the time before the thrust reached its steady 


state value. Again the results were excellent: the sum of 
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the squares difference was only 1.52 degrees Kelvin. A plot 
of the experimental and model temperatures is shown in Figure 


1s. 
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Figure 15 Case 3: Experimental and Model Temperatures Vs. 
Time 

The values obtained for the unknown variables were 
ay=-0.2000, a,=5.9382, ay=100.00, b,=1.6236, and b,=0.0500. 
The variables lead to resistance values of R,=11.7085, 
R,=19.0253, and Rg=-0.6380. Again, the negative resistance of 
R, indicating heating from ground. The values of the 
convection heat transfer coefficients were calculated to be 
4786.95 W/mK and 114.26 W/mK at the tip and fin respectively. 

The values for b,, and b, found from NODE418.FOR were 


added to the geometric and material data in COEF.FOR to again 
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find the convective heat transfer coefficients at every time 
step. At the tip, the heat transfer coefficient decreased 
from an initial value of 6451.38 to the final value of 
4786.95. The heat transfer coefficient for the fin also 
decreased, from an initial value of 154.11 to the final value 
of 114.26. There was a 26% decrease at both the tip and fin, 
with the tip showing nonlinearities. The coefficients are 


plotted in Figure 16. 


18% Al CASE 


TIME (sec} 


Figure 16 Convective Heat Transfer Coefficients Plotted vs. 
Time For Case 3. 
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D. Erosion Front Modeling 

An energy balance equation can be written between the 
leading edge erosion heat flux, g/A,, and the heat required to 
maintain the vane leading edge ablation rate, or 


C( Ty~Ty) 


25 9pUnCp,(Tay7Ty) = SP peF [1+ 
A, - F 


] (27) 


where S, is the Stanton number, Ty is the leading edge 
recovery temperature, Ty, is the vane leading edge temperature, 
T, is the melting temperature of the vane material, F is the 
heat of fusion for tungsten, and C is the heat capacity of 


tungsten. Also note that 


Soh .U.Ce_=Ape (28) 


where Hy, the leading edge convection heat transfer 
coefficient, is found by a parameter identification program 
like one of those previously described. [Ref. 10:p. 2,3] 

A theoretical erosion rate can be found by manipulating 
equations (27) and (28) 


Hz ( Tay~ Ty) 


C(T,-T,) (29) 
Pref {1+———— | 


S= 


Ty can be estimated by running a four node simulation model 


and using the node one temperatures at each time step. 
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Equation (27) is based upon ablation of the vane, which 
requires that Ty>Ty. Therefore the erosion rate was set equal 
to zero until T, reaches T,. The melting temperature for the 
vane, which is a 90% tungsten-10% copper alloy by weight, is 
3513K. This temperature is higher than T, for both the 9% and 
18% cases, and therefore theoretically the vane should not 
erode. Since the vane does erode, Ty for the vane was taken 
as the melting temperature of copper, 1358K. This seemed 
reasonable since the melting point of copper is lower than 
that of tungsten. 

Once the erosion rate is found, it can then be integrated 
over the time of the firing to find a theoretical length of 
the vane eroded. This was done in both the 9% and 18% cases 
and is shown in Figures 17 and 18. The length of the vane 
eroded using this method is estimated as 1.1 cm for the 9% 
case and 2.3 cm for the 18% case. Although the 1.1 cm found 
for the 9% case is high compared to the 0.4 cm found 
experimentally, the 2.3 cm found for the 18% case is very 
close to the 2.5 cm found experimentally. 

Equation (29) can also be used to try and validate the use 
of the melting temperature of copper for T,. This was done by 
plotting the vane temperatures found in the simulation 
programs as a function of the length between nodes one and 
two, then using the known total length eroded from the 
experiment to find the apparent melting temperature. The 


plots for the 9% and 18% cases are shown in Figures 19 and 20. 
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Erosion rate, (cm/s) 


TIME (sec) 


Total length eroded (cm) 


TIME (sec) 


Figure 17. Erosion Rate and Total Length Eroded for the 9% 


Case 
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Erosion rate, (cm/s) 


TIME (sec) 


Total length eroded (em) 
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Figure 18 Erosion Rate and Total Length Froded for the 18% 
Case 
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The melting temperature found for the 9% case was 1732K while 
the Mareing temperature found for the 18% case was 1580K. 
Although both of these are higher than the melting temperature 
of copper, they are fairly close. The reason for the melting 
temperature of the vane being higher than predicted is due to 


the presence of tungsten which has a melting temperature of 


3683K. 
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Figure 19 Temperature Profiles Between Nodes One and Two For 
the 9% Case 


33 


TIME (200) 


~_ 
a 
2 
» 
a 
he 

3 
» 
= 
we 
= 
~ 
bed 
ay 
Qu 
a 
bl 
io 


Figure 20 Temperature Profiles Between Nodes One and Two for 
the 18% Case 
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Iv. DISCUSSION OF RESULTS 

The full scale three node model attempted to show that the 
three fin nodes of the Parker five node full scale model could 
be reduced to one node. This was done, obtaining similar 
results for the convection heat transfer coefficients at the 
tip and fin. This validated the use of only one fin node in 
Reno’s four node quarter scale model. 

The erosion effects of aluminized propellent on the 
quarter scale vanes had to be investigated. There were three 
postulates considered of how the heat transfer coefficients 
changed: 

(1) the heat transfer coefficients were independent of 
erosion rate and time, 

(2) the heat transfer coefficients were dependent upon 
erosion rate, but given a fixed erosion rate, were time 
independent, and 

(3) the heat transfer coefficients were dependent upon 
erosion rate and were time varient. 

The first postulate was investigated by A. Danielson in 
(Ref. 2]. He found that as the percentage of aluminum in the 
exhaust and the erosion rate increased, nonlinear factors 
began to have a larger impact and show the limitations of the 


linear model [Ref. 2:p. 9]. 
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To investigate the remaining two postulates, a model for 
the erosion of the vanes had to be developed. The erosion of 
the vanes was modeled as a linear decrease of the geometric 
dimensions as a function of time and mass loss percentage. 
This was done in the subprogram COEF. 

For the second and third postulates, the coefficients in 
the PID subprogram COEF were set to the appropriate values for 
cases two and three, thereby allowing the geometric dimensions 
to vary. This led to excellent results which remained fairly 
constant even as the percentage of aluminum in the propellant 
increased. The sum of the squares error was only 1.19 for the 
0% Al case, 0.73 for the 9% Al case, and 1.52 for the 18% Al 
case. This seemed to link the erosion rate to the heat 
transfer coefficients. 

To determine whether the heat transfer coefficients were 
time dependent, the program COEF.FOR was modified to calculate 
the heat transfer coefficients as a function of time. 
Although the heat transfer coefficients remained fairly 
constant at the fin, they decreased over time at the tip. 

An equation based on ablation of the vane was used to try 
to predict the erosion rate. The erosion rate was then 
integrated over the time of the motor firing to obtain the 
theoretical length of the vane which eroded. Although the 9% 
case predicted an eroded length which was more than double the 
experimental value, the 18% case was very close. Two of the 


reasons the 9% case was off can be explained by the simplicity 


36 


of the model and the assumption that ablation would being 
occuring at the melting temperature of copper instead of the 
tungsten-copper alloy which the vane was composed of. 

To find a closer value to the melting temperature of the 
vane, the simulated vane temperature was plotted as a function 
of length between nodes one and two. By using the known 
length of vane eroded, a theoretical melting temperature could 
be found. The melting temperatures found were 1732K and 1580K 
for the 9% and 18% cases respectively. This was much closer 
to the 1358K for the melting temperature of copper than the 


3513K for the tungsten-copper alloy of the vane. 
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CONCLUSIONS 


The five node full scale model can be reduced to a three 
node full scale model by removing two of the three fin 
nodes and produce comparable convective heat transfer 
coefficients. 


Erosion of thrust vector control vanes can be adequately 
modeled by a linear decrease of the geometric properties 
as a function of time and the percentage of aluminum used 
in the propellant. 


The negative values found for R, indicate heating of the 
vane from the mount area. 


Both the tip and fin convective heat transfer coefficients 
were dependant upon erosion rate and were time variant. 


The erosion rate and therefore the length of the vane 
which will erode over the time of a motor firing can be 
adequately predicted using an energy balance equation 
based upon ablation of the vane. 


The melting temperature of ‘the vane appears to be much 


closer to that of copper than the tungsten-copper alloy 
which is expected. 
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RECOMMENDATIONS FOR FURTHER STUDY 


® The erosion front modeling needs to be investigated 


further to see if erosion mechanisms other than ablation 
can be modeled such as direct impingement of the 
aluminized particles on the vane. 


The G-law erosion algorithm explained in [Ref. 9] may 
provide a method to use results from a quarter scale model 
to predict full scale heat transfer characteristics. 


The quarter scale model needs to be modified to include 
the heating effects in the vane mount area. 
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APPENDIX A. SIMULATION PROGRAM 


This appendix contains the FORTRAN code used in the 
program SIM.FOR, which is a forward model program to simulate 
the temperatures of a three node full scale model, and 
SIM4.FOR which is a forward model program to simulate the 


temperatures of the four node quarter scale models. 
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Program SIM 


This is a forward model program to simulate the 
temperatures of a three node full scale model. 


integer maxparam, neg 
parameter (maxparam=50,neq=3) 
integer idOd, istep, nout 


real*8 t,tend,a(3,3),b(3,3),u(3),t2(61),y(3) 
real*8 param(maxparam) ,fcon,float,a3q 


intrinsic float 

external fen, divprk, sset 
common/datal/a,b,u 

Open files for data input/output 


open(9,name='sim3.mat’, status=’new’ ) 
open(8,name='datam.dat’, status='old’) 


read in experimental data 


do i=1,61 
read(8,*) t2(i) 

enddo 

close (8) 


initialize matrices 


~2936 
-0147 
.0107 
0243 


oooo 


Hou bon 
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a3g=0.0001 
b(1,1)=1.0000 
b(2,2)=0.0500 


a(1l,1)=-(a(1,2)+b(1,1)) 
a(2,2)=- (a (2, 1)+a(2, 3) +b(2,2)) 
a(3,3)=- (a(3,2)+a3q) 
u(1)=2670 

u(2)=#2570 

u(3)= 


set initial conditions 
t=0.0 
do i=1,3 
y(i)=0.0 
enddo 


tol=0.0005 
call sset(maxparam, 0.0, param, 1) 


idod=1 

do istep=1,61 
tend=0.0768*float (istep) 
call DIVPRK(id0,neq,fcn,t,tend,tol,param,y) 
write (9,9001) t,t2(istep),y 

enddo 

final call to release workspace 


id0=3 
call DIVPRK(id0,neg,fen,t,tend,tol,param, y) 


9001 format (1£6.3,4f10.4) 
close (9) 
end 
subroutine fcn(neq,t,y,yprime) 
integer neq 


real*8 t,y (neq) ,yprime (neq) 
real*8 a(3,3),b(3,3),u(3),d 


common/datal/a,b,u 
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thrust profile simulation as step input 


if (t.gt.0.2) then 
d=1.0 

else 
d=0.0 

end if 


do i=1,neq 

yprime(i)=0.0 
do j=1,neq 
yprime (1) =yprime(i)+a(i,j)*y(j)+b(i,j) *u(j)*d 
enddo 

enddo 


return 
end 
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Program SIM4 


This is a forward model program to simulate the 
temperatures of a four node quarter scale model. 


integer maxparam,neq 
parameter (maxparam=50,neq=4) 
integer idO, istep, nout 


real*8 t,tend,a(4,4),b(4,4) ,u(4),y(4) 
real*8 param(maxparam) ,fcn,float,a4g 


intrinsic float 

external fcon, divprk, sset, coef 
common/datal/a,b,u 

Open files for data input/output 
open(9,name='sim49.mat’, status=’new’ ) 


initialize matrices 


do i=1,4 
do j=1,4 
a(i,j)=0.0 
b(i,j)=0.0 
enddo 
enddo 
u(1)=2155 
u(2)=2061 
u(3)=0.0 
u(4)=0.0 


set initial conditions 


t=0.0 

do i=1,4 
y(i)=0.0 

enddo 


tol=0.0005 
call sset(maxparam, 0.0, param, 1) 
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idod=1 
do istep=1,41 
tend=0.05*float (istep) 
call coef (tend) 
call DIVPRK(id0O,neq,fcen,t,tend,tol,param, y) 
write(9,9001) t,y 
enddo 
final call to release workspace 


id0=3 
call DIVPRK(id0,neq,fcn,t,tend,tol,param, y) 


9001 format (1f6.3,4£10.4) 
close (9) 
end 
subroutine fcn(neq,t,y,yprime) 
integer neg 


real*8 t,y(neq) , yprime (neq) 
real*8 a(4,4),b(4,4),u(4),d 


common/datal/a,b,u 
thrust profile simulation as step input 


if (t.gt.0.7) then 
1 


do i=1,neq 

yprime(i)=0.0 
do j=1,neq 
yprime (i) =yprime (i)+a(i,j)*y(j)+b(i,j)*u(j)*d 
enddo 

enddo 


return 
end 


subroutine coef (tend) 
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real*8 vt,vf,vs,atf,afs,1tf,lfs,rho,cp,k,sf 
real*8 vt0O,vf0,vs0,atf0,afs0,1tf0,1fs0,a12,a21,a23,a32 
real*8 a(4,4),b(4,4),u(4),¢c1,¢c2,¢c3,r12,r23 


common/datal/a,b,u 


vt0=2.6 
vf0=52.0 
vs0=23.0 


atf0=5.9 
afs0=5.2 


1tf0=5.0 
lfs0=6.0 


rho=18310.0 
cp=146.0 
k=173.0 
sf=0.25 


vt=vt0-0.0*tend 
vf=vf0-0.0*tend 
vs=vs0-0.0*tend 


atf=atf0-0.0*tend 
afs=afs0-0.0*tend 


ltf=1tf0-0.0*tend 
lfs=lfs0-0.0*tend 


r12=100.0*1tf£/ (k*atf) 
r23=100.0*1lfs/ (k*afs) 
cl=rho*cp*vt*0.000001 
c2=rho*cp*v£*0.000001 
c3=rho*cp*vs*0.000001 


rl2=r12/sf 
r23=r23/sf 
cl=ci*sfi**3 
c2Z=c2*si**3 
c3=c3*st**3 


a(1,2)=1/ (ci*r12) 
a(2,1)=1/ (c2*r12) 
a(2,3)=1/ (c2*r23) 
a(3,2)=1/ (c3*r23) 


a(3,4)=0.5376 


a(4,3)=0.1528 
a4g=-0.1651 
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b(1,1)=7.6511 
b(2,2)=0.0722 


return 
end 
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APPENDIX B. COEFFICIENT PROGRAM 


This appendix contains the FORTRAN code used in the 
program COEF.FOR which calculated the effect of erosion on the 
known coefficients of the A matrix and the heat transfer 


coefficients. 
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program coef 

integer i 

real*§ vt,vf,vs,atf,afs,ltf,lfs,t,rho,cp,k,sf 

real*8 vt0,vf0,vs0,atf0,afs0,1tf0,lfs0,a12,a21,a23,a32 
real*8 asf0,asf,ast0O,ast,b11,b22,ht,hf 


intrinsic float 


open (10,name=’ coef18.mat’,status=' new’) 
open(11,name='htc18.mat’,status=’new’ ) 


vt0=2.6 
vf£0=52.0 
vs0=23.0 


atf0=5.9 
afs0=5.2 


ast0=4.35 
asf0=112.16 


1tf0=5.0 
lfs0=6.0 


rho=18310.0 
cp=146.0 
k=173.0 
sf=#0.25 


b11=1.6236 
b22=0.05 


do i=1,33 
t=0.05*float (i) 
vt=evt0-0.8125*t 
vf=vf0-16.25*t 
vs=vs0-7.1875*t 


atf=atfoO-0.0*t 
afs=afs0-0.0*t 


ast=ast0-0.90625*t 
asf=asf0-23.367*t 


ltf=1tf0-1.5625*t 
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lis=lfis0-1.875*t 


y12=100.0*1tf£/ (k*atf£) 
r23=100.0*1fs/ (k*afs) 
cl=rho*cp*vt*0.000001 
c2=rho*tcp*vi*0.000001 
c3=rho*cp*vs*t0.000001 


rl2=r12/sft 
Y23=r23/st 
cl=cl*sf**3 
c2=c2*sfi*t*3 
c3=c3*sf**3 
ast=ast*sf**2 
asf=asf*sf**2 


al2=1/ (c1*r12) 
a2i=1/ (c2*r12) 
a23=1/ (c2*r23) 
a32=1/ (c3*r23) 


rf1l=1/ (b11*c1) 
rf£2=1/ (b22*c2) 


ht=10000.0/ (rf1*ast) 
hf=10000.0/ (rf2*asf) 


9999 format (1f10.4,2f10.2) 
9998 format (5f10.5) 


write (10,9998)t,a12,a21,a23,a32 
write (11,9999)t,ht,hf 

end do 

close (10) 

close (11) 


end 
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APPENDIX C. PID PROGRAMS 


This appendix contains the PID programs for the three node 
full scale model (NODE3), and the four node quarter scale 
models for propellant with 0% Al (NODE40), 9% Al (NODE49) and 


18% Al (NODE418). 
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Program NODE3 


This program is the PID program for the three node vane 


model. 


external temp 


integer m,n,iparm(6),ibtype,ldfjac 


parameter (m=61,n=3,ldfjac=m) 


real*8 rparm(7),x(n),f(m),xjac(m,n),xg(n),ssq,ubl1,ub2 
real*8 xlb(n),xub(n),xscale(n),fscale(m),float,ht,hf 
real*8 a(3,3),b6(3,3),u(3),t2(61),ys(3,61) 

real*8 rho,k,cp,sf,c1,¢c2,¢c3,r12,r23,a3g 

real*8 vt,vf,vs,atf,afs,ast,asf,ltf,lfs 


Variables 
m 
n 
iparm 
ibtype 
ldfjac 
rparm 
x(n) 
f (m) 
xjac (m,n) 


xg (n) 

xlb (n) 
xub (n) 
xscale(n) 


fscale(m) 


ssq 
a(3,3) 
b(3,3) 
u(3) 

t2 (61) 
ys (3,61) 
rho 

k 

cp 

ve 

vt 

vs 

att 


okt dat t bop tt 


fo op oan 


tot t Bod ub a bo opt te tt 


number of functions 

number of variables 

list of parameters for DBCLSF setup 
type of bounds on variables 

leading dimension of fjac 

list of parameters for DBCLSF setup 
the pt where the function is evaluated 
the computed function at the point x 
matrix containing a finite difference 
approx Jacobian at the approx solution 
initial guess of x 

x lower bound 

x upper bound 

vector containing the scaling matrix for 
the variables 

vector containing the scaling matrix for 
the functions 

sum of the squares 

a matrix 

b matrix 

[TRL, TR2, 0] 

experimental temperatures 

calculated temperatures 

density 

conduction heat transfer coefficient 
specific heat 

volume of the tip 

volume of the fin 

volume of the shaft 

cross sectional area from tip to fin 
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a 


afs = cross sectional area from fin to shaft 

ast = gurface area of the tip 

asf = Surface area of the fin 

LE = length from tip to fin 

lfs = length from fin to shaft 

sft = Scale factor 

ub1 = stagnation temperature, TR1 

ub2 = free stream temperature, TR2 

ht = convection heat transfer coefficient at 
tip 

hf = convection heat transfer coefficient at 
fin 


intrinsic float 
common/datal/a,b,u,t2,ys 


Open files for data input/output 


open(10,name=’result.dat’, status=‘new’ ) 
open(9,name=‘temp.mat’, status=’new’ ) 
open(8,name=’datam.dat’, status=’old’) 
open(7,name=’/input.dat’, status=’old') 


read in experimental data 


do i=1,61 
read(8,*) t2(i) 

enddo 

close (8) 


read in input data 


read (7,*) 
read(7,*) 

read (7,*) 
read(7,*) 
read(7,*) rho,k,ep 
read (7,*) 

read (7,*) 
read(7,*) vt, vf, vs 
read (7,*) 
read(7,*) 
read(7,*) atf, afs 
read (7,*) 

read(7,*) 

read(7,*) ast, asf 
read (7,*) 

read (7,*) 
read(7,*) 1ltf, lfs 
read (7, *) 
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read (7,*) 
read(7,*) sf, ubl, ub2 
close (7) 


initial conditions 
full scale data 


r12=100.0*1tf/ (k*at£) 
r23=100.0*1fs/ (k*afs) 
cl=rho*cp*vt*0.000001 
c2=rho*cp*v£*0.000001 
c3=rho*cp*vs*0.000001 


scaled data 


r12=r12/sf 
Y23=r23/st 
cl=c1l*st**3 
C2=c2*st**3 
c3=c3*st**3 


initialize matrices to zero 


do i=1,3 
u(i)=0.0 

do j=1,3 

a (1,7) =0:..0 

b(i,j)=0.0 

enddo 
enddo 
nh Eley acter 
a(2,1)=1/(c2*r12) 
a(2,3)=1/(c2*r23) 
aot 2) =1/ (c3*r23) 
a3g=0.0 
61,2) =0.0 
b(2,2)=0.0 
a(1,1)=-(a(1,2)+b(1,1)) 
a(2,2)=- (a(2,1)+a(2,3)+b(2,2) ) 
a(3,3)=-(a(3,2)+a3qg) 
u(1)=ub1 
u (2) =ub2 
xg (1) =a3g 
xg (2)=b(1,1) 
xgG(3)=b:C2 ,2) 
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set up parameters for DBCLSF call 


do i=l1,n 
xscale(i)=1.0 
xlb(i)=0.0001 
xub (i) =100.0 
xg (i)=0.01 
x(i)=0.0 

end do 


do i=1,m 
fscale(i)=1.0 

end do 

ibtype=0 


call dbclsf (temp,m,n,xg,ibtype,xlb,xub,xscale,fscale, 
& iparm, rparm,x,f,xjac,ldfjac) 


calculate unknown resistances and convection heat transfer 


coefficients 
a3gq =x (1) 
b(1,1)=x(2) 
b(2,2)=x(3) 


cl=rho*cp*vt*0.000001 
c2=rho*cp*vf£*0.000001 
c3=rho*cp*vs*0.000001 


Cl=cl*sf**3 
C2=c2*gt**3 
c3=c3*gi**3 

rf1 =1/(b(1,1)*c1) 
rf2 =1/(b(2,2)*c2) 
r3g =1/ (a3g*c3) 


ht =10000.0/ (rf1*ast) 
hf =10000.0/ (rf2*asf) 


print and save results 


write(6,*) ’ a3g bil b22’ 
write(6,9000) x(1),x(2),x(3) 


9000 format (3£12.4) 
9003 format (2£12.4) 


write(10,*)’ a3g b(1,1) b(2,2)’ 
write(10,9000) x(1),x(2),x(3) 
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write(10,*) 

write (10,*) 

write(10,*)’ rfl ri2 r3g’ 
write(10,9000) rfl1,rf2,r3g 

write (10,*) 

write(10,*) 

write(10,*)' ht hf’ 
write(10,9003) ht,hft 


write the temp-time data for MATLAB analysis 
do i=1,61 
tt=0.0768*float (i) 
write (9,9001)tt,ys(2,i),t2(i) 
enddo 
9001 format (1£6.2,2f10.3) 


close (10) 
close (9) 


end 

Subroutine TEMP (m,n,x,f) 

This calculates the temperature-time history using the 
current parameters supplied py DBCLSF called from PID. It 
calculates an error function returned to DBCLSF based on 
the differences between predicted and observed temperature 
histories. 

integer maxparam, neq 

parameter (maxparam=50, neq=3) 

integer id0,istep,nout,m,n 

real*8 t,tend,y(3),tol,fcn,float,param(maxparam) , 

real*8 x(3),£(61),coef 

real*8 a(3,3),b(3,3),u(3),t2(61),ys(3,61) 

real*8 rho,k,cp,sf,cl,c2,c3,r12,r23,a3g 

real*8 vt,vf,vs,atf,afs,ast,asf,1ltf,lfs 

intrinsic float 

external feon,divprk,sset 


common/datal/a,b,u,t2,ys 


open (1i2,name=’ incoming.dat’,status=’new’) 
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a3g =x (1) 
b(1, me 
b(2,2)=x(3) 


write (6,8000)a3g,b(1,1),b(2,2) 


8000 format (3£12.4) 


a(1,1)=- (a(1,2)+b(1,1)) 
a(2,2)=-(a(2,1)+a(2,3)+b(2,2)) 
a(3,3)=- (a(3,2)+a3q) 


Set initial conditions 


t=0.0 
do i=1,neq 
y(i)=0.0 
do j=1,61 
ys(i,j)=0.0 
enddo 
enddo 
tol=0.0005 


call sset (maxparm, 0.0, param, 1) 


id0d=1 
do istep=1,61 
tend=0.0768*float (istep) 
CALL DIVPRK (id0O, neq, fen, t, tend, tol, param, y) 
do i=1,3 
ys(i,istep) =y (i) 
enddo 
enddo 


Final call to release workspace 


idd=3 
call divprk (id0,neqg,fcn,t,tend,tol,param,y) 


calculate error functions 


do i=1,61 
f (i) =ys(2,i)-t2(i) 
enddo 


print out rms error 
ssqr=0.0 

do i=1,61 
ssqr=ssqr+f (i) *f (i) 
enddo 

ssqr=ssqr/61 
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xer=sqrt (ssqr) 

write(6,*) xer 

write (12,*) xer 
return 

end 


subroutine fen(neq,t,y,yprime) 
integer neq 


real*8 t,y(neq) , yprime (neq) 
real*8 a(3,3),b(3,3),u(3),d,ys(3,61) 


common/datal/a,b,u,t2,ys 
thrust profile simulation as step input 


if (t.gt.0.2) then 
d=1.0 

else 
d=0.0 

end if 


do i=1,neq 

yprime (i) =0.0 
do j=1,neq 
yprime (i) =yprime (i) +a(i,j) *y(j)+b(i,j) *u(j) *d 
enddo 

enddo 


return 
end 
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Program NODE40 


This program is the PID program for the four node vane 
model with ablation from exhaust with 0% Al. 


external temp 


integer m,n,iparm(6),ibtype,ldfjac 


parameter (m=122,n=5,1ldfjac=m) 


real*8 rparm(7),x(n),f(m),xjac(m,n),xg(n),ssq,ubl,ub2 
real*8 xlb(n),xub(n),xscale(n),fscale(m),float,ht,hf 
real*8 a(4,4),b(4,4),u(4),t3(61),t4(61),ys(4,61) 
real*8 rho,k,cp,sf,c1l,c2,c3,c4,r12,r23,a4qg 

real*8 vt0,vf0,vs0,atf0,afs0,ast0,asf0,1tf0,lfs0 
real*8 vt,vf,vs,atf,afs,ast,asf,1ltf,lfs 


Variables 
m 
n 
iparm 
ibtype 
ldfjac 
rparm 
x(n) 
£ (m) 
xjac (m,n) 


xg (n) 
xlb(n) 
xub (n) 
xscale(n) 


fscale(m) 


ssq 

a (neq, neq) 
b (neq,neq) 
u (neq) 

t3 (61) 

t4 (61) 

ys (neq, 61) 
rho 

k 

cp 

vt 

vt 


number of functions 

number of variables 

list of parameters for DBCLSF setup 
type of bounds on variables 

leading dimension of fjac 

list of parameters for DBCLSF setup 
the pt where the function is evaluated 
the computed function at the point x 
matrix containing a finite difference 
approx Jacobian at the approx solution 
initial guess of x 

x lower bound 

x upper bound 

vector containing the scaling matrix for 
the variables 

vector containing the scaling matrix for 
the functions 

sum of the squares 

a matrix 

b matrix 

{TR1, TR2, 0, 0] 

experimental temperatures at node 3 
experimental temperatures at node 4 
calculated temperatures 

density 

conduction heat transfer coefficient 
specific heat 

volume of the tip 

volume of the fin 
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vs = volume of the shaft 

att = cross sectional area from tip to fin 

afs = cross sectional area from fin to shaft 

ast = surface area of the tip 

asf = surface area of the fin 

Lee = length from tip to fin 

lfs = length from fin to shaft 

sf = scale factor 

ubl = stagnation temperature, TR1 

ub2 = free stream temperature, TR2 

ht = convection heat transfer coefficient at 
tip 

hf = convection heat transfer coefficient at 
fin 


intrinsic float 


common/datai/a,b,u,t3,t4,ys 
common/data2/rho,k,cp,sf,c1,c2,c3 
common/data3/vt0,vf0,vs0,atf0,afs0,ast0,asf0,1tf0,1fs0 
common/data4/vt,vf,vs,atf,afs,ast,asf,ltf,lfs 


Open files for data input/output 


open(10,name=’result0O.dat’, status=‘new’ ) 
open(9,name=’temp0.mat’, status=’new’ ) 
open(8,name=‘datam0.dat’', status=’o0l1d’) 
open(7,name=’input.dat’, status=’old’) 


read in experimental data 


do i=1,61 
read(8,*) t3(i) 

enddo 

do i=1,61 
read(8,*) t4(i) 

enddo 

close (8) 


read in input data 


read (7,*) 

read (7,*) 

read (7, *) 

read (7, *) 

read(7,*) xho,k,cp 
read(7,*) 

read (7,*) 

vtO, vf£0, vs0 
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read(7,*) atf0, afso 
read (7,*) 

read (7,*) 

read(7,*) ast0O, asf0 
read(7,*) 

read (7, *) 

read(7,*) 1tf0, lfso 
read (7,*) 

read(7,*) 

read(7,*) sf, ubl, ub2 
close (7) 


initial conditions 


t=0 
tend=0 
call coef (x,tend) 


set up parameters for DBCLSF call 


do i=1,n 
xscale(i)=1.0 
i (i) =-0.2 
b(i) oe 0 
(1) 
x(i)= 
end do 


. 


do i=1,m 
fscale(i)=1.0 
end do 


ibtype=0 


call dbclsf(temp,m,n,xg,ibtype,xlb,xub,xscale,fscale, 
& iparm, rparm,x,f,xjac,ldfjac) 


calculate unknown resistances and convection heat transfer 
coefficients 


a(3,4)=x(1) 
a(4,3)=x(2) 
a4g =x (3) 
b(1,1)=x(4) 
b(2,2)=x(5) 
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ci=rho*cp*vt*0.000001 
c2=rho*cp*v£*0.000001 
c3=rhotcp*vs*0.000001 


Ccl=cl*sfi**3 
c2=c2*sf**3 
c3=c3*sf**3 


rfl =1/(b(1,1) *c1) 
rvf2 =1/(b(2,2) *c2) 
Y34 =1/(a(3,4)*c3) 
c4 =1/(a(4,3)*r34) 
r4g =1/(a4g*c4) 


ht =10000.0/ (rf1* (ast*sf**2) ) 
hf =10000.0/ (rf2* (asf*sf**2) ) 


print and save results 


write(6,*) ’ a34 a43 a4g bli b22' 
write(6,9000) x(1),x(2),x(3),x(4),x(5) 


9000 format (5£10.4) 
9003 format (2£11.4) 


write (10,*)/’ a34 a43 a4g b(1,1) 
b(2,2)’ 

write(10,9000) x(1),x(2),x(3),x(4),x(5) 

write (10,*) 

write (10, *) 

write (10,*)’ rfl rf2 r4g' 
write(10,9000) rf1,rf2,r4g 

write (10, *) 

write (10, *) 

write (10,*)’ ht hf’ 
write(10,9003) ht,hf 


write the temp-time data for MATLAB analysis 


do i=1,61 
tt=0.05*float (i) 
write (9,9001)tt,ys(3,1),ys(4,1),t3(i),t4(i) 
enddo 
9001 format (2x,1f6.4,4£10.4) 


close (10) 
close (9) 


end 
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Subroutine TEMP (m,n,x,f) 


This calculates the temperature-time history using the 
current parameters supplied by DBCLSF called from PID. It 
calculates an error function returned to DBCLSF based on 
the differences between predicted and observed temperature 
histories. 


integer maxparam, neq 

parameter (maxparam=50, neq=4) 

integer id0,istep,nout,m,n 

real*8 t,tend,y(4),tol,fcn,float,param(maxparam) , 
real*8 x(n),f(m),coef 

real*8 a(4,4),b(4,4) ,u(4),t3(61),t4(61),ys(4,61) 
real*8 rho,k,cp,sf,cl1,c2,c3,c4,r12,r23,a4g 
real*8 vt0O,vf0,vs0,atf0,afs0,ast0,asf0,1tf0,lfs0 
real*8 vt,vf,vs,atf,afs,ast,asf,1tf,lfs 
intrinsic float 

external fcn,divprk,sset,coef 
common/datal/a,b,u,t3,t4,ys 
common/data2/rho,k,cp,sf,c1,c2,c3 
common/data3/vt0,vf0,vs0,atf0,afs0,ast0,asf0,1tf0,1fs0o 
common/data4/vt,vf,vs,atf,afs,ast,asf,ltf,lfs 


open (12,name=’ incoming.dat’,status=‘new’ ) 


a(3,4)=x(1) 
a(4,3)=x(2) 
a4g =x (3) 
b(1,1) =x(4) 
b(2,2) =x(5) 


write (6,8000)a(3,4),a(4,3),a4g,b(1,1),b(2,2) 
8000 format (5£10.4) 

a(1,1)=-(a(1,2)+b(1,1)) 

a(2,2)=-(a(2,1)+a(2,3)+b(2,2)) 

a(3,3)=-(a(3,2)+a(3,4)) 

a(4,4)=- (a(4,3)+a4g) 

Set initial conditions 

t=0.0 

do i=1,neq 
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y(i)=0.0 
do j=1,61 
ys(i,j)=0.0 
enddo 
enddo 


tol=0.0005 
call sset (maxparm, 0.0, param, 1) 


idod=1 
do istep=1,61 
tend=0.05*float (istep) 
call coef (x,tend) 
CALL DIVPRK (idO, neq, fen, t, tend, tol, param, y) 
do i=1,4 
ys(i,istep) =y(i) 
enddo 
enddo 


Final call to release workspace 


id0=3 
call divprk (id0,neqg,fcn,t,tend,tol,param,y) 


calculate error functions 


do i=1,61 
£(i)=ys(3,1i)-t3 (i) 
£ (1+61) =ys (4,1) -t4(i) 
enddo 


print out rms error 
ssqr=0.0 

do i=1,m 
ssqr=ssqrtf (i) *f (i) 
enddo 

ssqr=ssqr/m 
xer=sqrt (ssqr) 
write(6,*) xer 
write (12,*) xer 
return 

end 


subroutine fen(neq,t,y,yprime) 
integer neq 


real*8 t,y (neq) ,yprime (neq) 
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real*8 a(4,4),b(4,4),u(4),d,ys (4,61) 


common/datal/a,b,u,t3,t4,ys 
common/data2/rho,k,cp,sf,cl,¢c2,c¢3 
common/data3/vt0,vf0,vs0,atf0,afs0,ast0,asf0,1tf0,1lfs0o 
common/data4/vt,vf,vs,atf,afs,ast,asf,ltf,lfs 


thrust profile simulation as step input 


if (t.gt.0.3) then 
d=1.0 

else 
d=0.0 

end if 


do i=1,neq 

yprime (i) =0.0 
do j=1,neq 
yprime (i) =yprime (i) +a(i,j)*y(4)+b(i,j) *u(4) *d 
enddo 

enddo 


return 
end 


subroutine coef (x,tend) 

integer i,j 

real*8 tend,x(5)} 

real*8 a(4,4),b(4,4),u(4) 

real*8 rho,k,cp,sf,cl,c2,¢c3,c4,r12,7r23,a4g 
real*8 vt0,vf0,vs0,atf0,afs0,ast0,asf0,1tf0,1lfso 
real*8 vt,vf,vs,atf,afs,ast,asf,ltf,lfs 
common/datal/a,b,u,t3,t4,ys 
common/data2/rho,k,cp,sf,c1,c2,c3 
common/data3/vt0,vf0,vs0,atf0,afs0,ast0,asf0,1tf0,lfso0 
common/data4/vt,vf,vs,atf,afs,ast,asf,1tf,lfs 
a,b matrix modification due to ablation effects 
full scale data 

vt=vt0-0.013*tend 

vf=vf0-0.26*tend 

vs=vs0-0.115*tend 


atf=atf0-0.0*tend 
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afa=afs0-0.0*tend 


ast=ast0-0.0145*tend 
asf=asf0-0.374*tend 


ltf=1t£0-0.025*tend 
lfs=lfs0-0.03*tend 


v12=100.0*1t£/ (k*atf) 
123=100.0*1fs/ (k*afs) 
cl=rho*cp*vt*0.000001 
c2=rho*cp*vi*d.000001 
c3=rho*cp*vs*0.000001 


scaled data 


r1l2=r12/sf 
r23=r23/sf 
Cl=cl*sf**3 
C2=c2*sft*3 
c3=c3*sti**3 


a(1,2)=1/(cl*r12) 
a(2,1)=#1/ (c2*r12) 
a(2,3)2#1/(c2*r23) 
a(3,2)=1/ (c3*r23) 


a(3,4)=x(1) 

a(4,3)=x(2) 

a4g =x (3) 

b(1,1) =x(4) 

b(2,2) =x(5) 
a(1,1)=-(a(1,2)+b(1,1)) 
a(2,2)=- (a(2,1)+a(2,3)+b(2,2)) 
a(3,3)=-(a(3,2)+a(3,4)) 
a(4,4)=- (a(4,3)+a4qg) 


return 
end 
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Program NODE49 


This program is 


the PID program for the four node vane 


model with erosion from exhaust with 9% Al. 


external temp 


integer m,n,iparm(6),ibtype,ldfjac 


parameter (m=82, 


real*8 rparm(7), 


n=5,ldfjac=m) 


x(n),£(m),xjac(m,n),xg(n),ssq,ubl,ub2 


real*8 xlb(n),xub(n),xscale(n),fscale(m),float,ht,hf 
real*8 a(4,4),b(4,4),u(4),t3 (41) ,t4(41) ,ys (4,41) 
real*8 rho,k,cp,sf,cl,c2,¢c3,c4,r12,7r23,a4g 

real*8 vt0,vf0,vs0,atf0,afs0,ast0,asf0,1tf0,1lfs0 


real*8 vt,vf,vs, 


Variables 
m 
n 
iparm 
ibtype 
ldfjac 
rparm 
x(n) 
£(m) 
xjac (m,n) 


hi tt tot to tot 


xg (n) 

x1b (n) 

xub (n) 

xscale(n) 
th 

fscale (m) 
th 

ssq 

a (neq,neq) 

b (neq, neq) 

u (neq) 

t3 (41) 

t4 (41) 

ys (neq, 41) 

rho 

k 

cp 

vt 

ve 

vs 

att 

afs 


fete eaae ne OU OH Gg y 


foto to ono 


atf,afs,ast,asf,ltf,lfs 


number of functions 

number of variables 

list of parameters for DBCLSF setup 
type of bounds on variables 

leading dimension of fjac 

list of parameters for DBCLSF setup 
the pt where the function is evaluated 
the computed function at the point x 
matrix containing a finite difference 
approx Jacobian at the approx solution 
initial guess of x 

x lower bound 

x upper bound 

vector containing the scaling matrix for 
variables 

vector containing the scaling matrix for 
functions 

sum of the squares 

a matrix 

b matrix 

(TRL, TR2, 0, 0] 

experimental temperatures at node 3 
experimental temperatures at node 4 
calculated temperatures 

density 

conduction heat transfer coefficient 
specific heat 

volume of the tip 

volume of the fin 

volume of the shaft 

cross sectional area from tip to fin 
cross sectional area from fin to shaft 
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ast = surface area of the tip 

asf = surface area of the fin 

itf = length from tip to fin 

lfs = length from fin to shaft 

sf = scale factor 

ubl = stagnation temperature, TR1 

ub2 = free stream temperature, TR2 

ht = convection heat transfer coefficient at 
tip 

hf = convection heat transfer coefficient at 
fin 


intrinsic float 


common/datal/a,b,u,t3,t4,ys 
common/data2/rho,k,cp,sf,c1,c2,c3 
common/data3/vt0,v£0,vs0,atf0,afs0,ast0,asf0,1tf£0,1fs0 
common/data4/vt,vf,vs,atf,afs,ast,asf,ltf,lfs 


Open files for data input/output 


open(10,name=’result9.dat’, status=‘’new’ ) 
open(9,name=’temp9.mat’, status=’new’) 
open(8,name=’datam9.dat’, status='old’) 
open(7,name=’input.dat’, status=’ol1d’) 


read in experimental data 


do i=1,41 
read(8,*) t3(i) 

enddo 

do i=1,41 
read(8,*) t4(i) 

enddo 

close (8) 


read in input data 


read (7,*) 

read(7,*) 

read (7,*) 

read (7,*) 

read(7,*) xho,k,cp 
read (7,*) 

read (7,*) 

read(7,*) vtO, vf0, vso 
read (7,*) 

read(7,*) 

read(7,*) atf0, afso 
read(7,*) 

read (7, *) 
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read(7,*) ast0, asf0 
read (7,*) 

read(7,*) 

read(7,*) 1tf0, lfso 
read(7,*) 

read(7,*) 

read(7,*) sf, ubl, ub2 
close (7) 


initial conditions 


t=0 
tend=0 
call coef (x, tend) 


u(1)=ubl 
u(2) =ub2 
u(3)=0.0 
u(4)=0.0 


set up parameters for DBCLSF call 


do i=1,n 
xscale(i)=1.0 
xlb(i)=-0.2 
xub (i) =100.0 
xg(i)=0.1 
x(i)=0.0 

end do 


do i=1,m 
fscale(i)=1.0 

end do 

ibtype=0 


call dbclisf(temp,m,n,xg,ibtype,xlb,xub,xscale,fscale, 
& iparm, rparm,x,f,xjac,ldfjac) 


calculate unknown resistances and convection heat transfer 
coefficients 


a(3,4)=x(1) 
a(4,3)=x(2) 
a4g =x (3) 


b(1,1) =x(4) 
b(2,2)=x(5) 


cl=rho*cp*vt*O.000001 


c2=rho*cp*v£*0.000001 
c3=rho*cp*vs*0.000001 
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Cl=ci*st**3 
e2=c2*st**3 
c3=c3*st**3 


rfl =1/ (b(1,1)*c1) 

rf2 =1/(b(2,2)*c2) 

034 =1/(a(3,4)*c3) 

c4 =1/(a(4,3)*r34) 
r4g =1/ (a4g*c4) 


ht =10000.0/ (rfi* (ast*sf**2) ) 
hf =10000.0/ (rf2* (asf*sf**2) ) 


print and save results 


write(6,*) ’ a34 a43 a4g bil b22!' 
write (6,9000) x(1),x(2),x(3),x(4),x(5) 


9000 format (5£10.4) 
9003 format (5x,2£10.4) 


write(10,*)’ a34 a43 a4g b(1,1) b(2,2)' 
write(10,9000) x(1),x(2),x(3),x(4),x(5) 

write (10, *) 

write (10,*) 

write(10,*)’ rfl rf2 rag’ 
write(10,9000) rf1,rf2,r4qg 

write (10, *) 

write (10, *) 

write (10,*)’ ht hf’ 

write(10,9003) ht,hf 


write the temp-time data for MATLAB analysis 
do i=1,41 
tt=0.05*float (i) 
write (9,9001)tt,ys(3,1),ys(4,1),t3(1),t4(1) 
enddo 
9001 format (2x,1£6.4,4£10.4) 


close(10) 
close (9) 


end 
Subroutine TEMP (m,n,x,f) 
This calculates the temperature-time history using the 


current parameters supplied by DBCLSF called from PID. It 
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calculates an error function returned to DBCLSF based on 
the differences between predicted and observed temperature 
histories. 


integer maxparam, neq 

parameter (maxparam=50, neq=4) 

integer id0,istep,nout,m,n 
real*8t,tend,y(4),tol,fcn,float,param(50),x(n),f£(m),coef 
real*8 a(4,4),b(4,4),u(4),t3(41),t4(41),ys(4,41) 
real*8 rho,k,cp,sf,cl,c2,c3,c4,r12,r23,a4g 

real*8 vt0,vf0,vs0,atf0,afs0,ast0O,asf0,1tf0,lfso 
real*8 vt,vf,vs,atf,afs,ast,asf,ltf,lfs 

intrinsic float 

external fcen,divprk,sset,coef 
common/datal/a,b,u,t3,t4,ys 
common/data2/rho,k,cp,sf,c1,c2,¢c3 
common/data3/vt0,vf0,vs0,atf0,afs0,ast0,asf0,1tf0,lfs0 
common/data4/vt,vfi,vs,atf,afs,ast,asf,ltf,lfs 


open (12,name=’ incoming9.dat’,status=’new’ ) 


write (6,8000)a(3,4),a(4,3),a4g,b(1,1) ,b(2,2) 


8000 format (5£10.4) 


a(1l,1)=-(a(1,2)+b(1,1)) 
a(2,2)=-(a(2,1)+a(2,3)+b(2,2)) 
a(3,3)=-(a(3,2)+a(3,4)) 
a(4,4)=- (a(4,3)+a4g) 


Set initial conditions 


t=0.0 
do i=1,neq 
y(i)=0.0 
do j=1,41 
ys(i,j)=0.0 
enddo 
enddo 
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tol=0.0005 
call sset (maxparm, 0.0, param, 1) 


idO=1 
do istep=1,41 
tend=0.05*float (istep) 
call coef (x, tend) 
CALL DIVPRK (idO, neq, fen, t, tend, tol, param, y) 
do i=1,4 
ys(i,istep) =y(i) 
enddo 
enddo 


Final call to release workspace 


id0=3 
call divprk (id0,neq,fcn,t,tend,tol,param,y) 


calculate error functions 


do i=1,41 

f (i) =ys(3,1i)-t3(i) 

£ (i+41) =ys (4,1) -t4(i) 
enddo 


print out rms error 

ssqr=0.0 

do i=i,m 
ssqr=ssqrtf (i) *f(i) 

enddo 

ssqr=ssqr/m 

xer=sqrt (ssqr) 

write(6,*) xer 

write(12,*) xer 

return 

end 


subroutine fcn(neq,t,y,yprime) 
integer neq 


real*8 t,y(neq) , yprime (neq) 
real*8 a(4,4),b(4,4),u(4),d,ys(4,41) 


common/datal/a,b,u,t3,t4,ys 
common/data2/rho,k,cp,sf,c1,c2,c3 
common/data3/vt0,vf0,vs0,atf0,afs0,ast0,asf0,1tf0,lfs0 
common/data4/vt,vf,vs,atf,afs,ast,asf,ltf,lfs 
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thrust profile simulation as step input 


if (t.gt.0.7) then 
d=1.0 

else 
d=0.0 

end if 


do i=1,neq 

yprime (i) =0.0 
do j=1,neq 
yprime (i)=yprime(i)+a(i,j)*y(j)+b(i,j) *u(j)*d 
enddo 

enddo 


return 
end 


subroutine coef (x, tend) 
integer i,j 


real*8 tend,x(5) 

real*8 a(4,4),b(4,4),u(4) 

real*8 rho,k,cp,sf,cl,c2,c3,c4,r12,r23,a4g 
real*8 vt0O,vf0,vs0,atf0,afs0,ast0,asf0,1tf0,1lfs0 
real*8 vt,vf,vs,atf,afs,ast,asf,1ltf,lfs 


common/datal/a,b,u,t3,t4,ys 
common/data2/rho,k,cp,sf,c1l,c2,c3 
common/data3/vt0,vf0,vs0,atf0,afs0,ast0,asf0,1tf0,lfs0 
common/data4/vt,vf,vs,atf,afs,ast,asf,1tf,lfs 

a,b matrix modification due to ablation effects 

full scale data 

vt=vt0-0.0*tend 

vf=vf0-0.0*tend 

vs=vs0-0.0*tend 


atf=atf0-0.0*tend 
afs=afs0-0.0*tend 


ast=ast0-0.0*tend 
asf=asf0-0.0*tend 


ltf=1tfoO-0.0*tend 
lfs=lfs0-0.0*tend 
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r12=100.0*1tf£/ (k*atf£) 
£23=100.0*1lfs/ (k*tafs) 
ecl=rho*cp*vt*0.000001 
c2=rho*cp*v£*0.000001 
c3=rho*cp*vs*0.000001 


scaled data 


r12=r1i2/sf 
r23=r23/sft 
cl=c1l*sf**3 
c2=c2¥*si**3 
c3=c3*si**3 


a(1,2)=1/ (c1*r1i2) 
a(2,1)=1/ (c2*r12) 
a(2,3)=1/ (c2*r23) 
a(3,2)=1/ (c3*r23) 


a(3,4)=x(1) 

a(4,3)=x(2) 

a4g =x (3) 

b(1,1) =x (4) 

b(2,2)=x(5) 

a(1,1)=- (a(1,2)+b(1,1)) 
a(2,2)=- (a(2,1)+a(2,3)+b(2,2)) 
a(3,3)=- (a(3,2)+a(3,4)) 
a(4,4)=-(a(4,3) +a4g) 
return 

end 
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Program NODE418 


This program is the PID program for the four node vane 
model with ablation from exhaust with 18% Al. 


external temp 
integer m,n,iparm(6),ibtype,ldfjac 
parameter (m=66,n=5,ldfjac=m) 


real*8 rparm(7),x(n),f(m),xjac(m,n),xq(n),ssq,ub1,ub2 
real*8 xlb(n),xub(n),xscale(n),fscale(m),float,ht,hf 
real*8 a(4,4),b(4,4) ,u(4),t3(33),t4(33),ys(4,33) 
real*8 rho,k,cp,sf,cl,c2,c3,c4,r12,r23,a4g 

real*8 vt0,vf0,vs0,atf0,afs0,ast0,asf0,1tf0,lfs0 
real*8 vt,vf,vs,atf,afs,ast,asf,ltf,lfs 


Variables 

m = number of functions 

n = number of variables 

iparm = list of parameters for DBCLSF setup 

ibtype = type of bounds on variables 

ldfjac = leading dimension of fjac 

rparm = list of parameters for DBCLSF setup 

x(n) = the pt where the function is evaluated 

£ (m) = the computed function at the point x 

xjac(m,n) = matrix containing a finite difference 

approx Jacobian at the approx solution 

xg (n) = initial guess of x 

x1b (n) = x lower bound 

xub (n) = xX upper bound 

xscale(n) = vector containing the scaling matrix for 
the variables 

fscale(m) = vector containing the scaling matrix for 
the functions 

ssq = sum of the squares 

a(neq,neq) = a matrix 

b(neq,neq) = b matrix 

u (neq) = [TR1, TR2, 0, 0] 

t3 (33) = experimental temperatures at node 3 

t4(33) = experimental temperatures at node 4 

ys (neq,33) = calculated temperatures 

rho = density 

k = conduction heat transfer coefficient 

cp = specific heat 

vt = volume of the tip 

vf = volume of the fin 

vs = volume of the shaft 

att = cross sectional area from tip to fin 
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afs = cross sectional area from fin to shaft 

ast = surface area of the tip 

asf = surface area of the fin 

ltt = length from tip to fin 

lfs = length from fin to shaft 

sf = scale factor 

ubl = stagnation temperature, TRI 

ub2 = free stream temperature, TR2 

ht = convection heat transfer coefficient at 
tip 

hf = convection heat transfer coefficient at 
fin 


intrinsic float 


common/datal/a,b,u,t3,t4,ys 
common/data2/rho,k,cp,sf,¢1,¢2,c3 
common/data3/vt0,vf0,vs0,atf0,afs0,ast0O,asf0,1tf0,lfso 
common/data4/vt,vf,vs,atf,afs,ast,asf,ltf,lfs 


Open files for data input/output 


open(10,name='result18.dat’, status='new’ ) 
open(9,name=’temp18.mat’, status=’new’) 
open(8,name=’datam181.dat’, status=’old’) 
open(7,name='input.dat’, status='old’) 


read in experimental data 


do i=1,33 
read(8,*) t3(i) 

enddo 

do i=1,33 
read(8,*) t4(i) 

enddo 

close (8) 


read in input data 


read (7,*) 

read (7,*) 

read (7, *) 

read (7, *) 

read(7,*) rho,k,ep 
read(7,*) 

read (7, *) 

read(7,*) vt0O, v£0, vso 
read (7,*) 

read (7,*) 

read(7,*) atfO, afso 
read (7,*) 
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ast0, asf0o 


1tfo, lfso 
+) 

pea *) sf, ubl, ub2 
close(7) 


initial conditions 


t=0 

tend=0 

call coef (x, tend) 
u(1)=ubl1 
u(2)=ub2 
u(3)=0.0 
u(4)=0.0 


set up parameters for DBCLSF call 


do i=1,n 
xscale(i)=1.0 
xlb(i)=-0.2 
xub (i) =100.0 
xg (i)=0.1 
x(i)=0.0 

end do 


do i=1,m 


fscale(i)=1.0 
end do 


ibtype=0 


call dbclsf(temp,m,n,xg,ibtype,xlb,xub,xscale,fscale, 
& iparm, rparm,x,f,xjac,ldfjac) 


calculate unknown resistances and convection heat transfer 
coefficients 


a(3,4)=x(1) 
a(4,3)=x(2) 
a4g =x (3) 
b(1,1)=x(4) 


b(2,2)=x(5) 


ci=rho*cp*vt*0.000001 
c2=rho*cp*vi*d.000001 
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c3=rho*cp*vs*0.000001 


cl=cl*sf**3 
Cc2=c2*sf**3 
c3=c3*si**3 


rfl =1/(b(1,1) *cl) 

r£2 =1/(b(2,2)*c2) 

r34 =1/(a(3,4)*c3) 

c&4 =1/(a(4,3)*r34) 
r4g =1/ (a4g*c4) 


ht =10000.0/ (rf1* (ast*sf**2) ) 
hf =10000.0/ (rf2* (ast*sf**2) ) 


print and save results 


write(6,*) ' a34 a43 a4g bil b22’ 
write(6,9000) x(1),x(2),x(3),x(4),x(5) 


9000 format (5f10.4) 
9003 format (2x,2f10.4) 


write(10,*)’ a34 a43 a4g b(1,1) b(2,2)’ 
write(10,9000) x(1),x(2),x(3),x(4),x(5) 

write (10, *) 

write (10, *) 

write (10,*)’ rfl rf2 r4g’ 
write(10,9000) rf1,rf2,r4g 

write (10,*) 

write (10,*) 

write (10,*)’ ht hf’ 

write(10,9003) ht,hf 


write the temp-time data for MATLAB analysis 
do i=1,33 
tt=0.05*float (i) 
write (9,9001)tt,ys(3,i),ys(4,i),t3(i),t4(i) 
enddo 
9001 format (2x,1f6.4,4£10.4) 


close(10) 
close (9) 


end 


Subroutine TEMP (m,n,x,f£) 
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This calculates the temperature-time history using the 
current parameters supplied by DBCLSF called from PID, It 
calculates an error function returned to DBCLSF based on 
the differences between predicted and observed temperature 
histories. 


integer maxparam, neq 

parameter (maxparam=50, neq=4) 

integer id0O,istep,nout,m,n 

real*8 t,tend,y(4),tol,fen,float,param(maxparam) , 
real*8 x(n),f£(m),coef 

real*8 a(4,4),b(4,4),u(4),t3(33),t4(33),ys(4,33) 
real*8 rho,k,cp,sf,cl,c2,c3,c4,r12,r23,a4g 
real*8 vt0,vf0,vs0,atf0,afs0,ast0,asf0,1ltf0,lfs0o 
real*8 vt,vf,vs,atf,afs,ast,asf,1ltf,lfs 
intrinsic float 

external fon,divprk,sset,coef 
common/datai/a,b,u,t3,t4,ys 
common/data2/rho,k,cp,sf,cl,c2,¢3 
common/data3/vt0,vf0,vs0,atf0,afs0,ast0,asf0,1tf0,1lfs0 
common/data4/vt,vf,vs,atf,afs,ast,asf,ltf,lfs 


open (12,name='incoming18.dat’, status=’ new’ ) 


a(3,4)=x(1) 
a(4,3)=x(2) 
a4g =x (3) 
b(1,1) =x(4) 
b(2,2)=x(5) 


write (6,8000)a(3,4),a(4,3),a4g,b(1,1),b(2,2) 
8000 format (5£10.4) 


-(a(1,2)+b(1,1)) 
-(a(2,1)+a(2,3)+b(2,2)) 
~({a(3,2)+a(3,4)) 
-(a(4,3)+a4g) 


pp @ 


(1,1) 
(2,2) 
(3,3) 
(4,4) 


, 
, 
, 
, 


Hoot 


Set initial conditions 


t=0.0 
do i=1,neq 
y(i)=0.0 

do j=1,33 
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ys(i,j)=#0.0 
enddo 
enddo 


tol=0.0005 
call sset (maxparm, 0.0, param, 1) 


id0=1 
do istep=1,33 
tend=0.05*float (istep) 
call coef (x,tend) 
CALL DIVPRK (idO, neq, fen, t, tend, tol, param, y) 
do i=1,4 
ys (i,istep) =y (i) 
enddo 
enddo 


Final call to release workspace 


idod=3 
call divprk (id0,neq,fcen,t,tend,tol,param,y) 


calculate error functions 


do i=1,33 

f£ (i) =ys (3,1) -t3 (i) 
£(1+33) =ys (4,1) -t4(i) 
enddo 


print out rms error 
ssqr=0.0 

do i=1l,m 
ssqr=ssqr+f (i) *£ (i) 
enddo 

ssqr=ssqr/m 
xer=sqrt (ssqr) 
write(6,*) xer 
write(12,*) xer 
return 

end 


subroutine fen(neq,t,y,yprime) 
integer neq 


real*8 t,y(neq) , yprime (neq) 
real*8 a(4,4),b(4,4),u(4),4,ys (4,33) 
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common/datai/a,b,u,t3,t4,ys 
common/data2/rho,k,cp,sf,cl,c2,c3 
common/data3/vt0,vf£0,vs0,atf0,afs0,ast0,asf0,1tf0,1fso 
common/data4/vt,vf,vs,atf,afs,ast,asf,ltf,lfs 


thrust profile simulation as step input 


if (t.gt.0.1) then 
d=1.0 

else 
d=0.0 

end if 


do i=1,neq 

yprime (i) =0.0 
do j=1,neq 
yprime (i) =yprime (i) +a(i,j)*y(j)+b(i,j) *u(j)*d 
enddo 

enddo 


return 
end 


subroutine coef (x, tend) 
integer i,j 


real*8 tend,x(5) 

real*8 a(4,4),b(4,4),u(4) 

real*8 rho,k,cp,sf,c1,c2,¢3,¢c4,r12,r23,a4g 
real*s8 vt0,vf0,vs0,atf0,afs0,ast0,asf0,1tf0,1lfs0 
real*8 vt,vf,vs,atf,afs,ast,asf,ltf,1fs 


common/datal/a,b,u,t3,t4,ys 
common/data2/rho,k,cp,sf,cl1,c2,c3 
common/data3/vt0,vf£0,vs0,atf0,afs0,ast0O,asf0,1tf0,lfso 
common/data4/vt,vfi,vs,atf,afs,ast,asf,1ltf,lfs 

a,b matrix modification due to ablation effects 


full scale data 
vt=vt0-0.0*tend 
vf=vf0-0.0*tend 
vs=vs0-0.0*tend 


atf=atrf0-0.0*tend 
afs=afs0-0.0*tend 
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ast=-ast0-0.0*tend 
asf=-asf0-0.0*tend 


ltf=ltf0-0.0*tend 
lfs=lfs0-0.0*tend 


y12=100.0*1tf£/ (k*atf) 
123=100.0*1lfs/ (k*afs) 
cl=rho*cp*vt*d0.000001 
c2=rho*cp*vf£*0.000001 
c3=rho*cp*vs*0.000001 


scaled data 


rl2=r12/st 
r23=r23/sf 
cl=cl*sf**3 
C2=c2*sir*3 
c3=c3*sf**3 


a(1l,2)=1/ (c1*r12) 
a(2,1)=1/ (c2*r12) 
a(2,3)=1/(c2*r23) 
a(3,2)=1/ (c3*r23) 


a(3,4)=x(1) 

a(4,3)=x(2) 

b(2,2)=x(5) 
a(1,1)=-(a(1,2)+b(1,1)) 
a(2,2)=- (a(2,1)+a(2,3)+b(2,2)) 
a(3,3)=-(a(3,2)+a(3,4)) 
a(4,4)=- (a(4,3)+a4q) 


return 
end 


82 


APPENDIX D. PHYSICAL DATA FILES 


The physical data files for the PID programs which 
contains the geometric and material properties of the vanes 


and the recovery temperatures used in each case. 
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NAWC INVERSE HEAT TRANSFER PROGRAM. INPUT DATA FOR NODE3.FOR 


Material properties: 


rho (kg/m*3), k (w/m-deg k), Cp (J/kg deg k) 
18310.0 173.0 146.0 
Vol (tip, cm*3), Vol (fin, cm‘*3), Vol (shaft, cm*3) 
2.6 52.00 23.0 
X-section areas: tip-fin (cm*2) fin-shaft (cm*2) 
539 5.2 
Surface areas: tip (cm*2) fin (cm*2) 
4.35 112.16 
Conductive path tip-fin (cm) fin-shaft (cm) 
5.0 6.0 
Scale factor: TR1 TR2 
1.00 2670 2570 
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NAWC INVERSE HEAT TRANSFER PROGRAM. INPUT DATA FOR NODE40.FOR 


Material properties: 


rho (kg/m‘*3), k (w/m-deg k), Cp (J/kg deg k) 
18310.0 173.0 146.0 
Vol (tip, cm*3), Vol (fin, cm*3), Vol (shaft, cm*3) 
2.6 52.00 23.0 
X-section areas: tip-fin (cm*2) fin-shaft (cm*2) 
5.9 BZ 
Surface areas: tip (cm*2) fin (cm*2) 
4.35 112.16 
Conductive path tip-fin (cm) fin-shaft (cm) 
5.0 6.0 
Scale factor: TRL TR2 
0.25 2360 2260 
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NAWC INVERSE HEAT TRANSFER PROGRAM. INPUT DATA FOR NODE49.FOR 


Material properties: 


rho (kg/m*3), k (w/m-deg k), Cp (d/kg deg k) 
18310.0 1735:0 146.0 
Vol (tip, cm*3), Vol (fin, cm*3), Vol (shaft, cm‘3) 
2.6 52.00 23.0 
X-section areas: tip-fin (cm*2) fin-shaft (cm*2) 
5.9 5.2 
Surface areas: tip (cm*2) fin (cm*2) 
4.35 112.16 
Conductive path tip-fin (cm) fin-shaft (cm) 
5.0 6.0 
Scale factor: TR1 TR2 
0.25 2464 2370 


86 


NAWC INVERSE HEAT TRANSFER PROGRAM. INPUT DATA FOR NODE418.FOR 


Material properties: 


rho (kg/m‘*3), k (w/m-deg k), Cp (J/kg deg k) 
18310.0 173.0 146.0 
Vol (tip, cm*3), Vol (fin, cm*3), Vol (shaft, cm*3) 
2.6 52.00 23.0 
X-section areas: tip-fin (cm*2) fin-shaft (cm*2) 
5.9 5.2 
Surface areas: tip (cm*2) fin (cm*2) 
4.35 112.16 
Conductive path tip-fin (cm) fin-shaft (cm) 
5.0 6.0 
Scale factor: TRL TR2 
0.25 2970 2870 
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APPENDIX E. IMSL ROUTINES 
A description of the IMSL routines DBCLSF, DIVPRK, and 


SSET used in the PID and simulation programs. 
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-BCLSF/DBCLSF (Single/Double precision} 


Purpose: Solve a nonlinear least squares problem subject <o bounds 
on the variables using a modified Levenberg-Marauardt 
algorithm and a finite-difference Jacobian. 

Usage: CALL BCLSF (FCN, M, N, XGUESS, IBTYPE, XLB, XUB, XSCALE, 

FSCALE, IPARAM, RPARAM, X, FVEC, FJAC, 
LDF JAC) 
Arguments 
FCN - User-supplied SUBROUTINE to evaluate the function to be 
minimized. The usage is 
‘ CALL FCN (M, N, -X, F), where 
M - Length of F. (Input) 
N - Length of X. (Input) 
x ~ The point at which the function is evaluated. 
(Input) 
X should not be changed by FCN. 
F - The computed function at the point X. 
(Output) 
FCN must be declared EXTERNAL in the calling program. 
M - Number of functions. (Input) 
N ~ Number of variables. (Input) 
XGUESS - Vector of length N containing the initial guess. (Input) 
TBTYPE ~- Scalar indicating the types of bounds on variables. 
(Input) 
IBTYPE Action 
Q User will supply all the bounds. 
Bi All variables are nonnegative. 
2 All variables are nonpositive. 
3 User supplies only the bounds on 1st variable, 
all other variables will have the same bounds. 
XLB - Vector of length N containing the lower bounds on 
variables. (Input, if IBTYPE = 0; output, if IBTYPE = 1 
or 2; input/output, if IBTYPE = 3) 
XUB - Vector of length N containing the upper bounds on 
variables. (Input, if IBTYPE = 0; output, if IBTYPE = 1 
or 2; input/output, if IBTYPE = 3) 
XSCALE - Vector of length N containing the diagonal scaling matrix 
for the variables. (Input) 
in the absence of other information, set all entries to 
1.0. 
FSCALE - Vector of length M containing the diagonal scaling matrix 
BCLSF/DBCLSF IMSL MATH/LIBRARY 
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for the functions. (Input) 
In the absence of other information, set ali entries to 
1.0. . 
IPARAM - Parameter vector of length 6. (Input/Output) 
See Remarks. 
RPARAM - Parameters vector of length 7. (Input/Output) 
See Remarks. 


x ~ Vector of length N containing the approximate solution. 
(Output) 
FVEC -~ Vector of length M containing the residuals at the 


approximate solution. (Output) 

FJAC - M by N matrix containing a finite difference approximate 
Jacobian at the approximate solution. (Output) 

LDFJAC - Leading dimension of FJAC exactly as specified in the 
dimension statement of the calling program. (Input) 


Remarks 


1. Automatic vorkspace usage is 
BCLSF 14*N + 2M - 1 units, or 
DBCLSF 26"N + 49M - 2 units. 
Workspace may be explicitly provided, if desired, by use of 
B2LSF/DB2LSF. The reference is 
CALL B2LSF (FCN, M, N, XGUESS, IBTYPE, XLB, XUB, XSCALE, 
FSCALE, IPARAM, RPARAM, X, FVEC, FJAC, 
LDFJAC, WK, IWK) 
The additional arguments are as follows: 
WK - Work vector of length 128N + 2eM - 1. WK contains 
the following information on output: 
The second N locations contain the last step taken. 
The third N locations contain the last Gauss-Newton step. 
The fourth N locations contain an estimate of the 
gradient at the solution. 
IWK ~- Work vector of length 2*N containing the 
permutations used in the QR factorization of the Jacobian 
at, the solution. 


2. Informational errors 
Type Code 

3 1 Both the actual and predicted relative reductions in the 

- function are less than or equal to the relative function 
convergence tolerance. 

4 2 The iterates appear to be converging to a noncritical 
. point. 

4 3 Maximum number of iterations exceeded. 


IMSL MATH/LIBRARY BCLSF/DBCLSF 
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4 4 Maximum number of function evaluations exceeded. 
4 5 Five consecutive steps have been taken with the maximum 


tep length. 


3. The first stopping criterion for BCLSF occurs when the norm of the 
function is less chan the absolute function tolerance. The second 
stopping criterion occurs when the norm of the scaled gradient 1s 
less than the given gradient tolerance. The third stopping 

riterion for BCLSF occurs when the scaled distance between the 
last two steps is less than the step tolerance. 


4. If nondefault parameters are desired for IPARAM or RPARAM, then 
U4LSF is called and the corresponding parameters are set to the 


desired value before 


calling the optimization program. Otherwise, 


if the default parameters are desired, then set IPARAM(1) to zero 
and call the optimization program omitting the call to U4LSF.. 
The call to U4LSF would be as follows: 


CALL U4LSF (¢ 


The following is a li 
IPARAM - Integer vect 
IPARAM(1) = 
IPARAM(2) = 


IPARAM(3) 
IPARAM(4) 
IPARAM(5) = 


IPARAM(6) = 


RPARAM - Real vector 
RPARAM(1) = 


RPARAM(2) = 
RPARAM(3) = 


RPARAM(4) = 


RPARAM(5) 
RPARAM(S) = 


BCLSF/DBCLSF 


IPARAM, RPARAM) . 


st of the parameters and the default values: 
or of length 6. 

Initialization flag. (0) 

Number of good digits in the function. 
(Machine dependent) 

Maximum number of iterations. (100) 

Maximum number of function evaluations. (400) 
Maximum number of Jacobian evaluations. (100) 
(Not used in BCLSF.) : 

Internal variable scaling flag. (1) 

If IPARAM(6) = 1 the values for XSCALE are 
set internally. 

of length 7. 

Scaled gradient tolerance. 

(SQRT(eps) in single precision) 

(eps**(1/3) in double precision) 

Scaled step tolerance. (eps*=(2/3)) 

Relative function tolerance. 
(MAX(1.0E-10,eps**(2/3)) in single precision) 
(MAX (1.0D-20,eps=*(2/3)) in double precision) 
Absolute function tolerance. 
(MAX(1.0E-20,eps=*2) in single precision) 
(MAX(1.0D-40,eps*=2) in double precision) 
False convergence tolerance. (100*eps) 
Maximum allowable step size. 
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= (1000*MAX(TOL1,TOL2)) where, 
TOL1 = SORT(sum of (XSCALE(I)=*XGUESS(I)) #2) 
for I= 1,...,N 
TOL2 = Z-norm of XSCALE. 
RPARAM(7) = Size of initial trust region radius. 
(Based on the initial scaled Cauchy step) 


eps is machine epsilon. 
If double precision is desired, then DU4LSF is called and RPARAM 


is declared double precision. 


Keyvords: Levenberg-Marquardt; Trust region 


Algorithm 
BCLSF uses a modified Levenberg-Marquardt method and an active set strategy to 
solve nonlinear jeast squares problems subject to simple bounds on the variables. 
The problem is stated as follows: 


td 1< 2 
main 5F(z)7 F(z) = 5 30 f(z)” 


wat) 
subject to /<r<u. 
where m > n. F : IR" — IR™. and f,(z) is the t-th component function of F(x). 
From a given starting point. an active set IA, which contains the indices of the 
variables at their bounds. is built. A variable is called a ‘free variable’ if it is not in 
the active set. The routine then computes the search direction for the free variables 
according to the formula 


d= —(J7I + wl) UTE, 


where y is the Levenberg-Marquardt parameter, F = F(z), and J is the Jacobian 
with respect to the free variables. The search direction for the variables in 1A is 
set to zero. The trust region approach discussed by Dennis and Schnabel (1983) 
is used to find the new point. Finally. the optimality conditions are checked. The 


conditions are: 
lolze)lse, Ua< uy 


92) <0, a= uy; 
(zi) > 0, ay = bh. 


where ¢ is a gradient tolerance. This process is repeated until the optimality criterion 
is achieved. 

The active set is changed only when a free variable hits its bounds during an 
iteration. or the optimality condition is met for the free variables but not for all 
_variables in IA. the active set. In the latter case, a variable which violates the 
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~ optimality condition will be dropped out of LA. For more detail on the Levenberg- 
Marquardt method. see Levenberg (1944). or Marquardt (1963). For more detailed 
information on active set strategy, see Gill and Murray {1976). 


Since a finite-difference method is used to estimate the Jacobian. for some single 
precision calculations. an inaccurate estimate of the Jacobian may cause the algo- 
rithm to terminate at a noncritical point. In such cases high precision arithmetic 
is recommended. Also. whenever the exact Jacobian can be easily provided. [MSL 
routine BCLSJ should be used instead. 


Example 
The nonlinear least squares problem 
RR 
min 5 x Atzy? 


com? 2 iat 


subject to ~2 


where f\(z) = 10(z2 — 27), and fo(z) = (1 — 21) is solved with an initial guess 
(—1.2, 1.0). and default values for parameters. 


¢ Declaration of variables 


INTEGER LDFJAC, M, N 
PARAMETER (LDFJAC#2, M=2, N=2) 


¢ 
INTEGER IPARAM(7), ITP, NOUT 
REAL FJAC(LDFJAC,N), FSCALE(M), FVEC(M), ROSBCK, 
& RPARAM(7), X(N), XGUESS(N), XLB(N), XS(N), XUB(N) 
EXTERNAL  BCLSF, ROSBCK, UMACH 
c Compute the least squares for the 
c Rosenbrock function. 
DATA XGUESS/-1.2E0, 1.0E0/, XS/2*1.0E0/, FSCALE/2#1.0E0/ 
DATA XLB/~2.0E0, -1.0£0/, XUB/O.5E0, 2.0E0/ 
c All the bounds are provided 
ITP * 0 ; 
c Default parameters are used 
IPARAM(1) # 0 
c 
CALL BCLSF (ROSBCK, M, N, XGUESS, ITP, XLB, XUB, XS, FSCALE, 
& IPARAM, RPARAM, X, FVEC, FJAC, LDFJAC) 
c Print results 
CALL UMACH (2, NOUT) 
WRITE (NOQUT,99999) X, FVEC, IPARAM(3), IPARAM(4) 
c 
99999 FORMAT (’ The solution 18 ’, 2F9.4, //, ' The function ’, 
& _@valuated at the solution is ‘, /, 18%, 2F9.4, //, 
& ’ The number of iterations is ', 10X, 73, /. ' The’, 
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& ‘number of function evaluations is ’, 13, /) 


END 
c 
SUBROUTINE ROSBCK (M, N, %, F) 
INTEGER M,N 
REAL X(N), FCM) 
c 
FC1) = 2. 0E1*(X(2)-X(1) =X(1)) 
F(2) = 1.0E0 - X(1) 
RETURN 
END 
Output 
The solution as - $000 - 2800 


The function evaluated at the solution is 
- 0000 «5000 


The number of iterations 1s 18 
The number of function evaluations is 22 
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- IVPRE/DIVPRK iSinale “Double precision | 


Purpose: Solve an initial-value problem for ordinary differential 
equations using the Runge-kutva-Verner fifth-order and 
Sixtn-order method. 
Usage. CALL IVPRK (IDO, NEQ, FCN. X, XEND, TOL, PARAN, Y) 
Arguments 
IDO - Flag indicating the state of the computation. 
(Input/Output) 
1 Initial entry 
2 Normal reentry 
3 Final call to release workspace 
4 Return’ because of interrupt 1 . 
5 Return because of interrupt 2 with step accepted 
6 _ Return because of interrupt 2 with step rejected 
Normally, the initial call is made with IDO=1. The 
routine then sets 1D0=2 and this value is then used for 
ail but the last call which is made with IDO=3. Thas 
final call is only used to release workspace, which was 
automatically allocated by the initial call with IDO=1. 
See Remark 2 for a description of the interrupts 
NEQ - Number of differential equations. (Input) 
FCK - User-supplied SUBROUTINE to evaluate functions. 
The usage is 
CALL FCN (NEQ, X, Y, YPRIME), where 
NEQ - Number of equations. (Input) 
x - Independent variable. (Input) 
Y - Array of length NEQ containing the dependent 
variable values. (Input) 
YPRIME - Array of length NEQ containing the vaiues of 
dY/aX at (X,Y). (Output) 
FCN must be declared EXTERNAL in the calling program 
. x ~ Independent variable. (Input/Output) 
On 2unput, X supplies the initial value. 
On output, X is replaced by XEND unless error conditions 
arise. See IDO for details. 
XEND - Value of X at which the solution is desired. (Input) 
XEND may be less than the initial value of X. 
TOL - Tolerance for error control. (Input) 
An attempt is made to control the norm of the local error 
such that the global error is proportional to TOL. 
More than one run, with different values of TOL, can be 
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usec to estimate the global error. 
Generally, 21% should not be greater than 0.001. 


' 


PARAM 


(Input/Output) 
If a parameter 1s zero then a default value is used. 
The Zollowing parameters must be set by the user. 


IVPRA.'DIVPRK 


PARAM 


ur 


Qo 


wo 


MXSTEP 


MXFCN 


INTRP1 


INTRP2 


- SCALE 


ITNORM 


Veczor of length SO containing optional parameters. 


Meaning 
Initial value of the step size H. 
Default: See Algorithm section. 
Minimum value of the step size H. 
Default: 0.0 
Maximum value of the step size KH. 
Default: No limit is imposed on the 
step size. 
Maximum number of steps allowed. 
Default: 500 
Maximum number of function evaluations 
allowed. 
Default: No limit 
Not used. 
£ nonzero then return with IDO=4. 
before each step. 
See Remark 3. 
Defauit: 0. 
If nonzero then return with IDO=5, 
after every successful step and with 
IDO=6 after every unsuccessful step. 
See Remark 3. 
Default: 0. 
A measure of the scale of the problem, 
such as an approximation to the average 
value of a norm of the Jacobian along 
the trajectory. 
Default: 1.0 
Switch determining error norm. 
In the following Ei is the absolute 
value of an estimate of the error in 
Y(i), called Yi here. 
O - min(absolute error. relative error) 
= max(Ei/Wi). i#1,2,...,NEQ, where 
Wi = max(abs(Yi.1.0), ms 
1 + absolute error = max(Ei), i=1,2... 
2 - max(Ei/Wi), i=1,2,..., where 
Wi = max(abs (Yi). FLOOR), 
and FLOOR is PARAM(11). 
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Remarks 


a? 


- Euclidian norm scaied ny YMAX 


= sartisum(E2**2/W142)). where 
wl = maxtabs(Yi),1.0). for YMAX. 
see Remark 3 
4: FLOOR - Usec in the norm computation 
Defauict: 1.0 
12-30 ~ Not used. 


The following entries in PARAM are set by the program. 


31 HTRIAL - Current trial step size 
32  4HMINC - Computed minimum step size allowed 
33 HMAXC - Computed maximum step size allowed. 
34 NSTEP - Number of steps taken. 
25 \NFCK - Number of function evaluations used. 
36-50 - Nov used. 

- Vector of length NEQ of dependent variables. 

(Input/Output) 


On input, Y contains the initial values. On output. 
Y contains the approxinate solution. 


i. Automatic workspace usage is 


IVPRK 
DIVPRK 


10=NEQ units. or 
20*NEQ units. 


Workspace may be explicitly provided, if desired, by use of 
The reference is 


I2PRK/DI2PRK 


CALL I2PRK (IDO, WEQ, FOX, X, XEND, TOL, PARAM, Y, 


VNORM, WK) 


The additional arguments are as follows: 
- User-supplied SUBROUTINE to compute the norm of the 


VNORM 


WK 


IMSL. Inc. 


error. 


(Input) 


The routine may be provided by the user, or the IMSL 
routine I3PRK/DI3PRK mer be used. 


The usage is 
CALL VNORM (NEQ, V, Y, YMAX, ENORM). where 


NEQ 
v 


bs 


YMAX 


ENORM 


Number of equations. (Input) 

Vector of length NEQ containing the vector whose 
norm 1s te be computed. (Input) 

Vector of length NEQ containing the values of 
the dependent variable. (Input) 

Vector of length NEQ containing the maximum Y 
values computed so far. (Input) 

Norm of the vector V. (Qutput) 


VNORM must be declared EXTERNAL in the calling program. 
- Work array of length i0*=NEQ. WK must not be changed 
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from the first call with IDO=1 until after the final call 


with IDO=3. 


Informational errors 

Type Code f 

4 1 Cannot satisfy error condition. TOL may be too small 
4 2 Too many function evaluations needed. 

4 3 Too many steps needed. The problem may be stiff. 


3. If PARAM(7) 1s nonzero. the subroutine returns with 
IDO = 4, and will resume calculation at the point of interruption 
if reentered with IDD = 4. If PARAM(8) is nonzero, the 
subroutine will interrupt the calculations 1mmediately after it 
decides whether or not to accept the result of the most 
recent trial step. IDO = 5 if the routine plans to accept, 
er IDD = 6 if it plans to reject. IDO may be changed by the user 
in order to force acceptance of a step (by changing IDO from 6 
to 5) that would otherwise be rejected. or vice versa. 
Relevant parameters to observe after return from an interrupt 
are IDO, HTRIAL, NSTEP, NFCN. and Y. Y is the newly computed 
trial value, accepted or not. 


Algorithm 


TVPRK ands an approximation to the solution of a system of first-order differential 
equations of the form y’ = f(z. y) with initial conditions. The routine attempts to 
keep the globai error proportional to a user-specified tolerance. The proportionality 
depends on the diderential equation and the range of integration. 

IVPRK is efficient for uonstitf systems where the derivative evaluations are not 
expensive and where the solution is not required at a large number of finely spaced 
points tas might be required for graphical output). 

IVPRK is based on a code designed by T. E. Hull, W. H. Enright and kK. R. 
Jackson /1976. 1977). It uses Runge-Kutta formulas of order five and six developed 
by J. H. Verner. : 


Example 


Consider a predator-prey problem with rabbits and foxes. Let r be the density 
of rabbits and let f be the density of foxes. In the absence of anv predator-prey 
interaction the rabbits would increase at a rate proportional to their number. and the 
foxes would die of starvation at a rate proportional] to their number. Mathematically. 


f= our 


fos =f. 
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The rate at which the rabbits are caten by une foxes is 2rf and the rate at whiel, 
the foxes increase, because they are cating the rabbits. is rf. So the mode} to be 


solved is 


ros 2r-Orf 
LU = -ferrf. 


The mitial conditions are r(Qi = 1 and f(O) = 3 over the interval 0 < 1 < 10. 

in the program Y(1) =r and Y(2) = f. Note that the parameter vector is firs: 
set to zero (using IMSL routine SSET} and then absolute error control is selected by 
setting PARAM(10) = 1.0. 

The last cal] to IVPRK with IDO = 3 releases IMSL workspace. that was reserved 
on the first call to IVPRK. It is not necessary to release the workspace in this 
example. because the program ends after solving a single problem. The call to 
release workspace is made as a model of what would be needed if the program 
included further calls to IMSL routines. 

The following piots are the result of using IVPRK with more closely spaced output 
than what is printed. (The program which does the plotting is not shown.) The 
second plot is a phase diagram for this system and clearly shows the periodic nature 
of the sojution. 


INTEGER MXPARM, NEQ 
PARAMETER (MXPARM=50, NEQ=2) 


c 
INTEGER IDO, ISTEP. NOUT 
REAL FCX, FLOAT. PARAM(MXPARM). T, TEND. TOL. Y(NEQ) 
INTRINSIC FLOAT 
EXTERNAL FCN. IVPRE, SSET, UMACH 
c 
CALL UMACH (2, NOUT) 
¢ Set anitial conditions 
T = 0.0 
YQ1) = 1.0 
¥(2) = 3.0 
c Set error tolerance 
TOL = 0.0005 
c Set -PARAM to default 
CALL SSET (MXPARM, 0.0, PARAM, 1) 
€ Select absolute error control 
PARAM(10) = 1.0 
c Print header 


WRITE (NOUT,99999) 
99999 FORMAT (4X, ‘ISTEP’. 5X, ‘Time’, 9X, ‘Y1', 11X. 'Y2’) 
IDO = 1 
DD 10 ISTEP#1. 10 
TEND = FLOAT(ISTEP) 
CALL IVPRK (IDO. NEQ. FCN, T, TEND. TOL, PARAM, Y) 
WRITE (NOUT, '(16,3F12.3)’) ISTEP, T, Y 
10 CONTINUE 
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c Final call to release workspace 


ID0 = 3 
“ CALL IVPRK (IDO. NEQ. FCN, T. TEND, TOL. PARAM. Y) 
END 


SUBROUTINE FCN (NEQ, T, ¥,. YPRIME) 
INTEGER NEQ 


REAL T. Y(NEQ), YPRIME(NEQ) 
c 
YPRIME(1) = 2,0¥(1) - 2.0=¥(1)*¥(2) 
YPRIME(2) = -¥(2) + Y(1)*¥(2) 
RETURN 
END 
Output 
ISTEP Time y1 Y2 
a§ 1.000 O78 1 465 
2 2.000 Oss 578 
3 3.000 . 292 250 
4 4.900 1.449 187 
5 5.000 4.046 1.444 
6 6.000 176 2.256 
7 7.000 . 066 308 
8 8.000 .148 367 
g 9.000 .655 188 
10 10.000 3.157 352 
4] 
A | Rabbits 
o } 
\ 
: \ ] 
ey \ | \ Foxes | 
| \ 
\ \ | 
\ 
i a 
KSZITSZ 
See Oe 
“9. 2. ‘. é. a. 10. 
Time 
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“Basic Linear Algebra Subprograms 


The basic linear algebra subprograms. usually called the BLAS, are routines for 
low-level vector operations such as dot products. Lawson et al. (1979) developed 
the original set of 38 BLAS routines. The IMSL BLAS collection includes these 
original 38 routines plus additional routines. The original BLAS are marked with 


a * in the descriptions. 


Programming Notes 


The BLAS do not follow the usual IMSL naming conventions. Instead the names 
consist of a prefix of one or more of the letters ‘I’. 'S'. ‘D’. °C’ and ‘2, a root 
name, and sometimes a suffix. For subprograms involving a mixture of data types 
the output type is indicated by the first prefix letter. The suffix denotes a variant 
algorithm. The prefix denotes the type of the operation according to the following 


table: 
I Integer 
S Real Cc 
D Double Z 
SD Single and double CZ 
DQ Double and quadruple ZQ 


Complex 

Double complex 

Single and double complex 
Double and quadruple complex 


Vector arguments have an increment parameter which specifies the storage space 
between elements. The correspondence between the vector z and the arguments SX 


and INC is 


SX((I-N) *INCX+1) 


if INCX < 0. 


as fe aeae es if INCX > 0 


Only positive values of INCX are allowed for operations which have a singie vector 


argument. 


The loops in all of the BLAS routines process the vector arguments in order of 
increasing i. For INCX < 0, this implies processing in reverse storage order. 


With the definitions, 


it 


Ui 


MX 
MY 
MZ 


max {1.1 +(N — 1)|INCX]} 
max {1.1 + (N = 1)/INCY]} 
= max {1,1 +(.V -1)[INcz|} 


the routine descriptions assume the following FORTRAN declarations: 


IMPLICIT INTEGER (I“N) 

IMPLICIT REAL s 

IMPLICIT DOUBLE PRECISION D 

IMPLICIT COMPLEX c . 
IMPLICIT DOUBLE COMPLEX 2 

INTEGER IX(MX) 

REAL SX(MX), SY(MY), SZ(MZ), SPARAM(S), 
& SH(LDH, =) 


BLAS 


102 


IMSL MATH/LIBRARY 


i 


Basic Linear Algebra Subprograms 


Integer | Real | Double Complex 


| Double 


| Complex | Page | 


' 
i 


rT, =a ;ISET  SSET | DSET CSET . ZSET 1026. 
y= fF, TTCOPY | SCOPY | DCOPY | CCOPY | ZCOPY 1026 | 
| m= ar, ,SSCAL | DSCAL | CSCAL  ‘ZSCAL ‘1026 | 
laéR | CSSCAL | ZDSCAL | 1026 | 
ly, = az, SVCAL | DVCAL | CVCAL | ZVCAL 1027 | 
'aélR | lesvcaL | zpvcaL | 1027 | 
[r=2,-0 ; ADD‘ SADD | DADD | CADD | ZADD | 1027, 
[n=a-ay | iSUB | SSUB | DSUB | CSUB 2SUB 1027 
Vi = az, +¥, 7 SAXPY | DAXPY | CAXPY | ZAXPY | 1027 
ry = 2, ISWAP | SSWAP | DSWAP | CSWAP ZSWAP | 1028 | 
ty | SDOT yopot” /cboTy | zpoTy | 1028 
| Fey | i !epote |2poTc ; 1028 
ir-yt ; DSDOT | CzbDOTU | zQDOTU | 1028 
Fey | czporc | zopotc | 1028 | 
ravr yt  SDSDOT | DODDOT | CZUDOT | ZOUDOT | 1028 
batty czcpoT | zocDoT | 1028 
;o+r-ut SDDOTI | DODOTI | CzboTI | zaboTI 1029 
‘acCeb—7-y 7 | sppoTa | DODOTA | CZDOTA | zapoTs | 1029 
Th = rey, | SHPROD | DHPROD | 1029 


1: 2, = min,z, 


SXYZ 


SPROCT 


IIMIN | ISMIN | IDMIN 


1: 7, max,2 TIMAX | ISMAX | IDMAX | 
| ICAMIN | IZAMIN 1032 


DXYZ 


DPRDCT 


rt: [ty, = max, [zj] TSAMAX | IDAMAX | ICAMAX | IZAMAX | 1031 
Construct Given’s | SROTG | DROTG | 1031 
rotation ! | | | 
Apply Given's | SROT DROT {| CSROT | 2DROT 1032 
rotation r 
Construct modified + SROTMG | DROTMG © 1032 

| Given's transform i | 
Apply modified SROTM | DROTM | CSROTM | ZDROTM | 1033 
Given’s transform 
Construct House- | SHOUTR | DHOUTR | 1034 
holder transform 
Apply Householder _ | SHOUAP | DHOUAP 1034 
transform - i i i 

tHigher precision accumulation used. 

BLAS 
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DOUBLE PRECISION DXCMX), DYC(MY), OZ2(MZ), DPARAM(5), 


& DH(LDH, «) 

DOUBLE PRECISION Dacc(2), DZACC (4) 
COMPLEX CXCMX), CYCMY) 
DOUBLE COMPLEX Z2X(MX), ZYCMY) 


Since FORTRAN 77 does not include the type DOUBLE COMPLEX. routines with 
DOUBLE COMPLEX arguments are not available for all systems. Some systems use 
the declaration COMPLEX*16 instead of DOUBLE COMPLEX. 

The set of BLAS routines are summarized by the table on page 1025. Routines 
marked with a dagger {}) in the table use higher precision accumulation. 


Set a Vector to a Constant Value 


CALL ISET (N, TA, IX, INCX) 
CALL SSET (N, SA, SX, INCX) 
CALL DSET (N, DA, DX, INCX) 
CALL CSET (N, CA, CX, INCX) 
CALL ZSET (N, ZA, 2X, INCX) 


These subroutines set x; = a for i = 1.2......N. [If V < 0 then the routines 
return immediately. 


Copy a Vector 
CALL ICOPY (N, [X, INCX, LY, INCY) 


“ CaLL SCOPY (N, SX, INCX, SY, INCY) 
* CALL DCOPY (N, DX, INCX, DY, INCY) 
* CALL CCOPY (N, CX, INCX, CY, INCY) 


CALL ZCOPY (N, ZX, INCX, ZY, INCY) 


These subroutines set yj = x; fori = 1.2.....M. If NV < 0 then the routines 
return immediately. 


Scale a Vector 


a CALL SSCAL (N, SA, SX, INCX) 
. CALL OSCAL (N, DA, DX, INCX) 
* CALL CSCAL (N, CA, CX, INCX) 
CALL ZSCAL (N, ZA, ZX, INCX) 
. CALL CSSCAL (N, SA, CX, INCX) 


CALL ZDSCAL (N, DA, 2X, INCX) 


These ‘subroutines set r; = az, fori = 1,2.....N. If N <$ 0 then the routines 
return immediately. 
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